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I.  INTRODUCTION 

This  thesis  addresses  efforts  to  improve  simulation,  using  conventional 
high  explosives  and  scale  models,  of  the  underwater  shock  environment  and 
structural  effects  resulting  from  the  underwater  nuclear  explosion. 
Simulation  using  scale  models  and  small  conventional  charges  provides 
valuable  information  without  the  requirement  for  nuclear  testing,  with 
minimal  environmental  impact,  and  at  low  cost.  Better  understanding  of  the 
physics  involved  impacts  ship  and  weapons  designs. 

In  the  current  atmosphere  of  reduced  military  spending  in  order  to  reap 
the  benefits  of  the  "Peace  Dividend"  resulting  from  the  end  of  the  cold  war, 
the  threat  to  U.  S.  Navy  ships  and  submarines  from  underwater  nuclear 
explosion  would  appear  to  be  greatly  reduced.  However,  two  factors  ensure 
the  continued  existence  of  the  threat  from  underwater  nuclear  explosion: 

1.  The  presence  of  the  former  Soviet  Union's  vast  arsenal  of  nuclear 
weapons  combines  with  economic  instability  to  increase  the  likelihood  of 
more  nations  gaining  access  to  the  material  necessary  to  construct 
nuclear  weapons. 

2.  The  ceaseless  march  of  technology  worldwide  dictates  future 
growth  in  the  number  of  nations  attaining  the  particular  technology 
necessary  to  build  and  detonate  nuclear  weapons. 

The  growing  community  of  nuclear  weapons  capable  nations  may  not  posses 

the  same  restraint  from  the  use  of  nuclear  weapons  displayed  since  1945. 

Given  the  existence  of  a  threat  from  underwater  nuclear  explosion,  the 

significance  of  this  threat  can  be  determined  only  if  the  effects  are  well 

understood.    This  same  understanding  is  essential  to  incorporation  of  shock 

hardening  in  ship  and  submarine  designs  to  improve  survivability. 


With  the  overall  objective  of  improving  ship  and  submarine 
survivability  through  better  understanding  of  the  phenomena  associated 
with  underwater  explosions,  the  Naval  Postgraduate  School  conducts 
ongoing  research  into  underwater  explosions  and  effects.  This  thesis  is  the 
result  of  a  part  of  that  continuing  research,  the  first  at  the  school  related 
specifically  to  the  nuclear  underwater  explosion. 

The  known  characteristics  of  the  nuclear  underwater  explosion, 
together  with  a  discussion  of  modeling  techniques  is  found  in  Chapter  II. 

Chapter  III  presents  a  comparison  of  the  structural  response  of  a  simple 
cylinder  subjected  to  side-on  attack  by  conventional  spherical  charge  and 
nuclear  type  pressure  profiles.  The  doubly  asymptotic  boundary  assumption 
combined  with  an  explicit  finite  element  method  was  used  to  perform  the 
analysis.  Results  verify  the  need  to  match  peak  pressure  and  duration  of  the 
nuclear  pressure-time  history  when  designing  test  charges  to  simulate  the 
underwater  nuclear  explosion. 

The  tapered  charge,  due  to  inherent  long  pressure  duration,  commonly 
generates  the  simulated  nuclear  pressure  field  used  in  model  testing. 
Chapter  IV  explores  the  three  dimensional  aspects  of  the  shock  front 
generated  by  a  conventional  underwater  tapered  charge  explosion.  The 
results  of  analysis  performed  using  an  explicit  finite  element  method  are 
presented  with  the  intent  of  supplementing  existing  knowledge  of  the 
tapered  charge  pressure  profile  which  to  date  consists  mostly  of  on-axis  data. 
The  results  show  the  evolving  shape  of  the  shock  front  at  early  times  and 
provide  information  on  the  relationship  between  the  peak  pressures 
measured  on  and  off  the  charge  axis. 


As  outlined  in  Chapter  V,  a  computer  code  was  written  to  optimize  the 
design  process  for  a  conventional  tapered  high  explosive  charge.  Starting 
with  a  desired  pressure-time  history  and  initial  estimates  for  charge 
geometry  and  standoff  distance,  this  program  utilizes  public  domain 
optimization  software  to  return  improved  design  values.  Although  tested 
with  a  subroutine  based  upon  an  existing  routine  for  calculating  the 
pressure-time  history  of  a  tapered  charge,  this  program  may  be  coupled  with 
other  existing  or  future  routines  for  calculating  the  tapered  charge 
pressure-time  history. 

Conclusions  and  recommendations  for  further  research  in  this  area  may 
be  found  in  Chapter  VI. 


II.  CHARACTERISTICS  AND  MODELING  OF  THE 
UNDERWATER  NUCLEAR  EXPLOSION 

Before  attempting  to  simulate  the  underwater  nuclear  explosion  using 
conventional  high  explosive  charges  for  scale  model  testing,  a  degree  of 
familiarity  with  the  shock  generated  by  nuclear  and  spherical  high  explosive 
charges  is  warranted.  This  chapter,  therefore,  outlines  the  use  of  empirical 
relationships  to  determine  the  pressure-time  histories  generated  by  the  two 
types  of  charges.  The  tremendous  weight  and  standoff  distance  required  for 
a  conventional  spherical  charge  to  create  a  pressure  profile  similar  to  that  of 
a  nuclear  charge  points  directly  to  the  need  for  scaling. 

Following  a  brief  description  of  the  principles  used  to  construct  scale 
models  for  underwater  shock  testing,  these  principles  are  applied  to  nuclear 
and  conventional  spherical  charges.  The  resultant  large  size  and  standoff, 
even  after  scaling,  of  the  conventional  spherical  charge,  leads  to  the  selection 
of  the  tapered  charge,  made  of  conventional  high  explosives,  to  simulate  the 
underwater  shock  generated  by  the  underwater  nuclear  explosion. 

A.  THE  UNDERWATER  NUCLEAR  SHOCK  PROFILE 

The  energy  content,  or  "yield",  of  a  nuclear  explosion  is  commonly 
measured  in  tons  of  TNT  equivalent,  the  amount  of  explosive  energy 
contained  in  2,000  pounds  of  the  conventional  high  explosive  TNT.  A  one 
kiloton  (kT)  device  contains  the  energy  equivalent  of  1,000  tons  of  TNT,  one 
megaton  (MT)  that  of  1,000,000  tons.  Although  the  majority  of  the  energy 
released  in  an  underwater  nuclear  explosion  contributes  to  the  generation  of 


the  underwater  shock  wave,  the  extremely  high  temperatures  (tens  of 
millions  of  degrees)  reached  in  a  nuclear  explosion  contribute  to  a  significant 
amount  of  energy  release  in  the  form  of  thermal  radiation.  Chemical 
explosions,  by  contrast,  occur  at  much  lower  temperatures  (thousands  of 
degrees),  resulting  in  a  higher  percentage  of  the  total  energy  released  as 
kinetic  energy  to  generate  the  underwater  shock.  (Glasstone  and  Dolan, 
1977,  pp.  1-3,6,  11) 

To  date,  the  United  States  has  conducted  five  announced  underwater 
nuclear  explosions,  from  1946  to  1962  (Bolt,  1976.  pp.  251-274).  Glasstone 
and  Dolan  (1977,  pp.  268-272)  provide  three  empirical  charts  to  calculate  the 
pressure-time  history  of  an  underwater  nuclear  explosion  given  the  yield  of 
the  device  and  the  standoff  distance  R  from  the  explosion.  From  the  curves 
of  the  first  two  charts,  the  maximum  pressure  Pmax  (psi)  and  the  time 
constant  8  (ms)  are  determined.  The  time  constant  equals,  as  in  exponential 
decay,  the  time  between  the  arrival  of  the  shock  when  P  =  Pmax  and  the  time 
at  which  P  =  Pmax/e  *  0.37Pmax.  The  pressure  actually  decays  at  a  somewhat 
slower  than  exponential  rate  after  one  time  constant,  necessitating  the  third 
chart  which  plots  the  non-dimensional  values  P(t)/Pmax  vs  t/0  to  provide  an 
idealized  pressure-time  history  for  an  underwater  nuclear  explosion  with  no 
bottom  or  surface  reflection  effects.  Based  upon  these  curves,  Figure  1  shows 
the  pressure  time  history  of  a  40  kT  nuclear  explosion  at  a  standoff  of  1,000 
yds.  Figure  2  illustrates  standoff.  Two  prominent  features  of  the 
underwater  nuclear  explosion  stand  out  in  Figure  1:  the  high  pressure  at  a 
significant  distance  from  the  explosion,   and  the  long  decay  time.     Pmax 


increases  with  yield,  decreases  with  standoff;  the  duration  of  the  shock  wave 
increases  with  yield  and  standoff  (Glasstone  and  Dolan,  1977,  p.  269). 
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Figure  1.  Pressure  profile  of  a  40  kT  nuclear  charge,  E  -  1000  yds. 
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Figure  2.  Measurement  of  standoff  distance  R. 
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In  water  of  sufficient  depth  where  bottom  reflections  are  negligible  or 
occur  at  a  much  later  time,  the  steady  pressure  decay  ends  abruptly  by  the 
phenomenon  of  surface  cutoff  (Shin  and  Geers,  1991,  §3.3;  or  Glasstone  and 
Dolan,  1977,  pp.  244-246).  Figure  3  shows  two  paths  followed  by  a  shock 
wave  emanating  from  an  underwater  explosion.  The  direct,  compression 
wave  strikes  the  target  first  with  a  sudden  rise  in  pressure  followed  by  the 
steady  pressure  decay  described  above.   The  other  path  shows  a  compression 


wave  striking  the  air- water  interface  and  being  reflected  as  a  tension,  or 
rarefaction  wave. 
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Figure  3.  Paths  followed  by  the  direct  and  rarefaction  waves. 

Approximating  shock  speed  by  the  speed  of  sound  in  water  Cs  gives  the 
cutoff  time  tc  as  the  difference  in  distance  traveled  by  the  direct  and 
rarefaction  waves  divided  by  C  : 


too  — 


2.    *   "  +  D2- R 


The  arrival  of  the  rarefaction  wave  at  time  tc  after  the  arrival  of  the 
direct  shock  wave  causes  a  sudden  pressure  drop  as  the  remaining  pressure 
from  the  compression  wave  is  essentially  canceled.     Hence,  if  the  40  kT 


nuclear  charge  discussed  earlier  is  detonated  at  a  285  ft  depth  in  deep  water, 
a  target  at  1,000  yards  and  the  same  depth  will  experience  t  =  10.7  ms  using 
the  expression  on  the  previous  page  with  Cs  =  5  ft/ms.  The  resultant 
pressure  profile  with  pressure  cutoff  is  shown  in  Figure  4. 
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Figure  4.  Pressure  profile  of  a  40  kT  nuclear  charge  with  surface  cutoff. 
R=1,000  yds,  D-285  ft. 

For  an  actual  underwater  nuclear  explosion,  the  smooth  curve  of  Figure 
4  would  be  further  modified  by  refraction  of  the  shock  wave  due  to  salinity, 
currents,  and  temperature  variation  in  the  water  media.  Whereas  Figure  4 
depicts  an  instantaneous  pressure  rise  upon  arrival  of  the  compression  wave 
at  time  0,  a  finite  rise  time  would  be  expected.  Additionally,  the  abrupt 
surface  cutoff  shown  assumes  the  rarefaction  wave  and  the  compression 
wave  travel  at  the  same  speed  of  sound  in  water.  In  actuality,  the 
rarefaction  wave  travelling  in  shocked  media  partially  overtakes  the 
compression  wave  rendering  a  less  steep  pressure  drop  at  tr( .  For  explosions 
occurring  in  shallow  water,  bottom  reflections  and  retransmissions  further 
alter  the  shock  profile.  Additional  shocks  may  occur  at  later  times  due  to  the 
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bubble  pulse  resulting  from  explosions  at  depths  such  that  the  gases  of  the 
explosion  expand  and  collapse  before  venting  at  the  surface.  More  than  three 
bubble  pulses  are  unlikely  due  to  steam  condensation.  (Glasstone  and  Dolan, 
1977,  pp.  56,  245,  246,  269) 

Having  established  the  general  characteristics  of  the  pressure-time 
history  of  the  underwater  nuclear  explosion,  the  next  section  uses  empirical 
formulas  to  determine  the  feasibility  of  emulating  this  profile  using 
conventional  high  explosives  of  spherical  shape. 

B.     SPHERICAL  CONVENTIONAL  CHARGE  PROFILE 

As  in  the  case  of  the  nuclear  charge,  the  conventional  charge  of 
spherical  or  near  spherical  shape  gives  rise  to  a  sudden  pressure  increase 
followed  by  exponential  pressure  decay  for  one  time  constant  9  and  somewhat 
slower  decay  thereafter.  Using  exponential  pressure  decay  for  an 
approximation,  empirical  studies  provide  the  following  useful  relationships 
for  determining  the  pressure-time  history  developed  by  a  conventional  high 
explosive  of  spherical  or  near  spherical  shape  when  detonated  underwater 
(Shin  and  Geers,  1991,  §3.2): 

P(t)  =  PmaXe:^ 

Pmax  =  Ki  11      =  maximum  pressure  (psi) 

0  =  K2W1/3(^-) ' '  =  time  (ms)  for  P  to  decay  to  %^ 

where  Ki,  K2,  Ai,  A2  =  empirical  constants  for  a  given  explosive 

t  =  time  (ms) 
P  =  pressure  (psi) 
R  =  standoff  (ft) 
W  =  charge  weight  (lb) . 
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Using  the  empirical  relationships  above,  a  charge  of  27.5  thousand  tons 
of  TNT  would  be  required  to  emulate  the  pressure-  time  history  generated  by 
a  40  kT  underwater  nuclear  explosion  as  illustrated  in  Figure  5. 
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Figure  5.   Comparison  of  40  kT  nuclear  (solid  line)  and  27,500  ton  TNT 
(dashed  line)  explosions.  R=3,000  ft,  D-285  ft  for  both  charges. 

The  fact  that  27.5  kT  of  TNT  generates  a  similar  pressure  profile  as  a 
nuclear  device  of  40  kT  yield  is  attributed  to  the  thermal  energy  released  by 
the  nuclear  explosion  as  discussed  in  Section  II.A.  Obviously,  using  27.5 
thousand  tons  of  TNT  to  study  the  effects  of  a  nuclear  explosion  is  not 
practical,  nor  feasible.  This  amount  of  high  explosive  represents  a  volume 
greater  than  that  of  the  Washington  Monument  (Lapp,  1980).  Scaling  laws, 
discussed  in  the  next  section,  enable  the  use  of  scale  models  and  charges  of  a 
more  reasonable  size  to  simulate  the  underwater  nuclear  explosion  and  its 
effects. 
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C.     HOPKINSON  SCALING 

Dimensional  analysis  yields  scaling  principles  which  enable  the  use  of 
scale  models  to  replicate  the  behavior  of  full  size  objects,  or  prototypes.  Of 
particular  utility  in  the  analysis  of  underwater  explosions  and  effects  is 
Hopkinson  scaling.  Through  the  use  a  scale  factor  X,  the  quantities  length, 
mass,  and  time  are  scaled  as  follows: 

Model  Prototype 

XL  Length  L 

Xi  Time  t 

X3M  Mass  M 

giving  the  invariant  quantities: 

p  Density  p 

c,  P  Stress,  Pressure  a,  P 

e  Strain  e. 

Some  quantities,  chiefly  the  hydrostatic  loading  due  to  gravity,   are  not 

adaptable  to  scaling  and  require  additional  consideration  in  modeling.   (Shin 

andGeers,  1991,  §4.1) 

The  benefits  of  scaling  are  many.  By  reducing  standoff  and  charge  size, 
tests  can  be  conducted  in  small  manmade  ponds  with  little  or  no 
environmental  impact.  Geometrically  similar  models  can  be  constructed  of 
the  same  materials  as  the  prototype  to  study  structural  response  to 
underwater  shock.  Model  stress  and  strain  levels  will  match  those  of  the 
prototype.  Since  the  material  required  to  build  model  charges  and  structures 
is  equal  to  X3  that  of  the  prototype,  model  testing  delivers  obvious  cost 
benefits. 
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Using  a  scale  factor  of  X  =  1/30  to  simulate  the  40  kT  nuclear 
pressure-time  history  discussed  previously  would  require  a  spherical  TNT 
charge  of  size 


WM  =  ?t3Wp  = 


30 


x28.5x  10  !  xtons  x2,000 


lbm 
ton 


2, 040  lbm. 


The  standoff  required  is 


RM  =  XRF  =  ^7*3, 000  ft  =  100  ft. 
Figure  6  shows  the  40  kT  nuclear  pressure-time  history  scaled  using  X  =  1/30 
together  with   the  pressure  profile  of  2,040  pounds   of  TNT.      The  only 
difference  between  this  figure  and  Figure  5  is  the  time  scale. 
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Figure  6.   2,040  lbm  TNT  charge  (dashed  line),  R=100  ft,  to  simulate  a 
40  kT  nuclear  (solid  line)  explosion  using  a  scale  factor  of  X-  1/30. 

Although  the  spherical  conventional  charge  size  is  now  feasible,  the 
weight  and  standoff  required  for  the  simulation  are  too  great  for  small  pond 
testing  limited  to  nominal  charge  weights  and  standoffs  in  the  tens  of  pounds 
and  feet  respectively.    One  might  be  tempted  in  scale  model  testing  to  use  a 
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smaller  charge  and  shorter  standoff  to  match  only  the  peak  pressure  of  the 
nuclear  profile  when  studying  shock  effects  on  a  scale  model.  As  shown  in 
Figure  7,  a  14  pound  spherical  TNT  charge  at  a  standoff  of  19  ft  gives  the 
same  maximum  pressure  as  a  2040  pound  charge  at  a  standoff  of  100  ft.  The 
pressure  of  the  lighter  charge  decays  much  more  rapidly  due  to  its  smaller 
time  constant. 
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Figure  7.  Pressure-time  histories  of  a  2,040  Ibm  (solid  line)  TNT  charge 
at  R=100  ft  with  a  14  Ibm  (dashed  line)  TNT  charge  at  R=19  ft. 

It  will  be  shown  in  Chapter  III  that  the  structural  response  from  two 
pressure  profiles  having  the  same  peak  pressure  but  different  duration  of 
high  pressure  give  rise  to  dramatically  different  structural  responses. 
Therefore,  in  order  to  simulate  the  slow  pressure  decay  of  the  underwater 
nuclear  explosion  using  conventional  charges  of  modest  size,  the  shape  of  the 
conventional  charge  must  be  modified.  The  resulting  shape  is  that  of  the 
tapered  charge  discussed  in  the  next  section. 
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D.  TAPED  CHARGE  PRESSURE  PROFILE 

Designed  to  generate  long  duration  shock  waves  to  simulate  nuclear 
underwater  shock  loading  on  scale  models,  tapered  charges  consist  of  a  series 
of  truncated  cones  on  a  common  axis  fitted  with  a  detonator  on  the  small  or 
nose  end  as  shown  in  Figure  8.  Constructed  in  sizes  ranging  from  a  few 
ounces  to  over  15,000  pounds,  the  tapered  charge  generates  a  directional 
pressure  field  with  maximum  duration  along  the  nose  side  on  the  charge 
axis.  (Gordon  and  Davidson,  1983) 


detonator 


Figxire  8.  Geometry  of  the  tapered  charge. 

For  tapered  and  spherical  charges  of  the  same  weight  and  at  the  same 
standoff,  Figure  9  illustrates  the  trade-off  between  peak  pressure  and 
duration  distinguishing  the  two  designs.  The  tapered  charge  creates  a  peak 
pressure  at  a  lower  value  than  the  spherical  charge,  followed  by  a  more 
gradual  pressure  decline  over  a  region  called  the  pressure  plateau.    At  the 
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end  of  the  pressure  plateau,  the  pressure-time  history  of  the  tapered  charge 
gives  way  to  exponential  decay  then  surface  cutoff. 
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Figure  9.  Pressure  profiles  of  tapered  (solid  line)  and  spherical  (dashed 
line)  charges  with  the  same  weight  and  standoff 

The  pressure  plateau  generated  by  the  tapered  charge  gives  an  obvious 
advantage  over  the  spherical  charge  in  simulating  the  nuclear  shock  profile. 
The  duration  of  the  pressure  plateau  t  can  be  determined  to  a  first 
approximation  by  estimating  the  time  difference  between  (1)  travel  from  nose 
to  tail  along  the  charge  axis  a  distance  L  at  the  detonation  speed  CD  then  a 
distance  L  +  R  to  the  target  at  the  speed  of  sound  in  water  Cs  ,  and  (2)  travel 
from  the  nose  to  the  target  a  distance  R  at  speed  Cs  .  The  resulting 
expression  for  pressure  plateau  duration  is 

R       r(  1  \\ 


tp  = 


L    |  L  +  R 

_Cd       Cs   j 


Cs=HcD  +  CsJ 
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The  above  formula  predicts  a  plateau  duration  independent  of  standoff.  As 
Figure  10  shows,  however,  the  duration  of  the  pressure  plateau  actually 
decreases  somewhat  with  increasing  standoff. 
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Figure  10.   Tapered  charge  pressure  profiles  at  varying  standoff. 

Figure  10  also  shows  the  approximate  1/R  relationship  between  peak 
pressure  and  standoff.  In  addition  to  varying  the  standoff  to  match  peak 
pressures  and  the  overall  length  to  match  pressure  plateau  duration,  the 
segment  lengths  and  diameters  of  the  tapered  charge  can  be  varied  to  taylor 
the  shape  of  the  tapered  charge  pressure-time  history  to  model  a  particular 
nuclear  pressure  profile.  Design  of  a  tapered  charge  involves  first  examining 
existing  data  for  a  charge  design  which  produces  a  pressure-time  history 
most  nearly  matching  that  desired.  The  charge  design  is  then  adjusted  using 
computer  calculations  tempered  with  the  experience  of  the  designer.    Small 
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scale  experiments  may  be  used  to  verify  the  design  before  production  of  a 
large  tapered  charge.  (Costanzo,  1991) 

Additional  knowledge  of  the  tapered  explosion  as  well  as  any 
streamlining  of  the  design  process  should  improve  the  design  of  tapered 
charges  used  to  simulate  the  underwater  nuclear  explosion.  An  application 
of  the  finite  element  method  (FEM)  to  study  the  early  time  propagation  of  the 
shock  generated  by  a  tapered  charge  explosion  is  found  in  Chapter  IV. 
Chapter  V  outlines  the  application  of  computer  design  optimization 
techniques  to  enhance  the  tapered  charge  design  process. 
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III.  EFFECT  OF  PRESSURE  DURATION  ON 
STRUCTURAL  RESPONSE 

As  mentioned  in  Section  II. C,  if  one  could  match  only  the  peak  pressure 
when  conducting  a  simulation  of  the  underwater  nuclear  attack,  a  spherical 
charge  of  modest  size  would  suffice.  As  illustrated  in  this  chapter,  however, 
matching  of  peak  pressures  alone  is  not  adequate.  Comparison  of  the 
response  of  a  simple  cylindrical  shell  to  side-on  attack  by  a  long  duration, 
nuclear  type,  and  a  short  duration,  conventional  type  pressure  profile  yielded 
substantially  different  results  using  numerical  techniques. 

A.     ATTACK  CURVES 

In  order  for  a  realistic  comparison  in  the  model  testing  environment,  the 
pressure  profiles  used  in  this  study  were  generated  from  two  56  lbm  HBX-1 
charges.  The  short  duration  pressure  profile  was  derived  from  the  empirical 
relationships  of  Section  II. B  for  a  spherical  charge  at  a  standoff  of  20  ft.  The 
long  duration  pressure  profile  used  in  this  study  to  simulate  the  nuclear 
profile  was  derived  from  scaled  tapered  charge  data.  Figure  11  shows  the 
two  attack  curves  with  the  same  peak  pressure  of  approximately  3,400  psi 
but  very  different  high  pressure  duration  times.  Each  pressure-time  history 
with  corresponding  standoff  was  entered  into  an  existing  computer  code  to 
provide  underwater  shock  loading  of  a  simple  cylindrical  aluminum  shell 
model  as  described  in  the  next  section. 
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Figure    11.      Two  attack   curves   used  for  comparison   of  structural 
response. 


B.     NUMERICAL  MODEL 

The  attack  geometry  for  this  study  is  shown  in  Figure  12.  The  figure 
illustrates  the  standoff  for  a  spherical  charge.  Tapered  charge  standoff  is,  by 
convention,  measured  from  the  nose,  Figure  8,  to  the  target.  The  pressure 
profiles  and  standoffs  from  the  previous  section  were  used  with  a  side-on 
attack  cylindrical  shell  model  developed  by  Fox  (1992,  pp.  60,61).  The 
analysis  was  conducted  using  the  public  domain  finite  element  method 
(FEM)  code  VEC/DYNA3D  (Hallquist  and  Stillman,  June,  1990)  coupled  with 
the  boundary  element  method  code  USA  (Deruntz,  1989).  The  USA 
(Underwater  Shock  Analyzer)  code  reduces  the  media  surrounding  the 
cylinder  and  the  associated  forces  to  discrete  forces  and  masses  to  provide 
loading  to  the  cylinder.  The  FEM  code  VEC/DYNA3D  utilizes  explicit  time 
integration  to  provide  the  response  of  the  cylinder  to  the  applied  shock 
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loading.      The   coupling   of  the    two    codes    was   initiated   by    the    Naval 
Postgraduate  School  (Fox,  Kwon,  and  Shin,  1992). 


R 
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Figure  12.  Attack  geometry  for  cylindrical  shell  FEM  experiment. 

The  target  consisted  of  a  1/4"  thick  6061-T6  aluminum  cylinder  with  1" 
machined  and  welded  endplates.  Material  properties  for  the  cylinder  were: 

oY  =  40,000  psi 

E  =  10,800  kpsi 

v  =  0.33 

p=  1741bm/ft3. 

Figure  13  depicts  the  550  shell  elements  comprising  the  one-quarter 
symmetry  discretization  of  the  cylindrical  shell  model  used  in  this  study.  A 
description  of  the  theory  behind  the  shell  elements  used  can  be  found  in  the 
article  by  Belytschko,  Lin,  and  Tsay  (1984).  Symmetric  boundary  conditions 
were  applied  on  the  yz  and  zx  planes.  The  boundary  element  loading 
described  previously  provided  the  boundary  conditions  for  the  outer  surface 
of  the  wet  elements  of  the  aluminum  cylinder. 
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Figure  13.  Quarter  symmetry  FEM  model  of  simple  cylindrical  shell. 

The  aluminum  was  modeled  as  a  kinematic/isotropic/elastic/plastic.  In 
this  idealization,  elastic  deformation  occurs  at  stress  levels  lower  than  the 
yield  stress  with  plastic  deformation  occurring:  at  the  yield  stress.  This 
provides  a  reasonable  approximation  of  the  uniaxial  stress-strain  curve  for 
6061-T6  aluminum  based  upon  pull-test  data  as  shown  in  Meyers  and  Murr 
(1981,  p.  40),  with  the  exception  of  no  provision  for  failure  once  a  maximum 
engineering  strain  of  7  to  9  per  cent  has  been  sustained.  Taking  8  per  cent 
plastic  strain  as  a  failure  criterion  and  applying  a  factor  of  safety  of  2,  failure 
of  the  aluminum  shell  model  was  predicted  for  effective  plastic  strain  in 
excess  of  4  per  cent    Effective  plastic  strain  ep  is  defined  as  (Ugural  and 

Fenster,  1987): 

J2 


eD  = 
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where  e,,  e2,  and  e3  are  the  true  plastic  strain  components.  Since  the 
maximum  strain-rate  calculated  in  this  study  was  approximately  100 
in/in-sec,  well  below  the  2000  in/in-sec  required  for  appreciable  strain-rate 
effects,  (Meyers  and  Murr.  1981,  p.  50),  strain-rate  hardening  was  not 
included  in  the  model. 

C.     RESULTS 

Figures  14  and  15  show  the  damage  sustained  by  the  cylindrical  target 
as  a  result  of  attack  by  the  nuclear  and  conventional  pressure  profiles  with 
the  same  peak  pressures.  The  quarter-symmetry  model  has  been  reflected 
about  the  yz  and  zx  planes  using  the  post-processor  TAURUS  (Hallquist, 
1990)  to  provide  visualization  of  the  full  model  in  each  case.  Displacements 
for  the  conventional  attack  of  Figure  14  were  scaled  by  a  factor  of  10  to  better 
display  the  damage  pattern.  No  scaling  was  done  in  Figure  15  where  the 
damage  was  much  more  extensive  due  to  the  nuclear  type  attack. 
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Figure  14.    Effective  plastic  strain  from  conventional  attack  (displace- 
ments scaled  by  a  factor  of  10). 
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Figure  15.   Effective  plastic  strain  from  nuclear  attack.    (No  scaling  of 
displacements) 


Figure  16  illustrates  the  time  history  of  effective  plastic  strain  for  the 
elements  sustaining  maximum  damage  from  each  attack.  The  maximum 
damage  of  1.67  per  cent  effective  plastic  strain  experienced  by  the 
conventional  attack,  although  significant,  did  not  exceed  the  failure  criteria 
selected  in  the  previous  section.  The  effective  plastic  strain  from  the  nuclear 
type  attack,  exceeding  4  per  cent  after  less  than  0.6  ms,  indicates  the 
prediction  of  catastrophic  failure  of  the  cylinder  resulting  from  this  attack 
based  upon  the  numerical  analysis. 

Not  only  is  the  magnitude  of  effective  plastic  strain  much  greater  for  the 
nuclear  type  attack,  the  mode  of  damage  is  quite  different.  In  both  cases, 
maximum  damage  occurred  in  elements  located  approximately  3.2  inches 
from  the  endplate.  In  the  conventional  case,  the  maximum  damage  occurred 
in  element  211,  on  the  yz  plane.  Maximum  damage  in  the  nuclear  case, 
however,  occurred  at  element  214,  30  degrees  off  the  yz  plane. 
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however,  occurred  at  element  214,  30  degrees  off  the  yz  plane.   The  locations 
of  the  two  critical  elements  are  shown  on  Figures  13  through  15. 


c 

(0 


10 
9 
8 
7 
6 
5 
4 
3 
2 
1 


* 



■-■■■  El  214,  Nuclear 

# 
* 

# 

i- 

1211. 

b 

Conventional 

• 

* 

1                • 

„.»' 

* 

_  ~ 

**^"^    *■ 

v0  '  1 



0.2 


0.4 


0.6 


0.8 


t(ms) 


Figure  16.    Plot  of  effective  plastic  strain  for  the  elements  sustaining 
maximum  damage. 


Based  upon  the  analysis  of  this  chapter,  both  the  peak  pressure  and 
high  pressure  duration  must  be  generated  by  the  conventional  high  explosive 
used  to  simulate  the  underwater  nuclear  explosion.  Due  to  the  ability  of  the 
tapered  charge  to  accomplish  this  task  at  a  lower  charge  weight  than  the 
spherical  charge,  the  studies  of  Chapters  IV  and  V  were  conducted  to 
enhance  understanding  and  improve  the  tapered  charge  design  process. 
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IV.  THE  TAPERED  CHARGE  SHOCK  FRONT 

Due  to  its  ability  to  develop  a  long  duration  pressure  profile  from  a 
modest  charge  weight,  the  tapered  charge  commonly  provides  the  shock 
loading  to  simulate  the  underwater  nuclear  attack  as  discussed  previously. 
In  order  to  better  understand  the  shock  developed  by  tapered  charges,  finite 
element  analysis  (FEA)  was  conducted  to  study  the  shock  developed  by  a 
tapered  charge  detonated  underwater.  The  FEA  provided  information  on  the 
directional  nature  of  the  pressure-time  history  generated  by  the  tapered 
charge.  This  information,  specific  for  the  charge  geometry  and  type  of 
explosive,  can  be  used  to  determine  the  accuracy  required  to  position  the 
charge  and  target  used  in  model  studies  of  the  underwater  nuclear  explosion. 

A.     FEA  MODEL 

Figure  17  shows  the  46"  x  46"  x  152"  quarter-symmetry  FEA  model  used 
for  this  study.  The  mesh  consists  of  120  HBX-1  charge  elements  and  40.692 
water  elements.  Appendix  A  provides  a  listing  of  the  input  file  for  the 
IN  GRID  (Stillman  and  Hallquist,  1991)  mesh  generating  program  used.  The 
dimensions  of  the  tapered  charge,  corresponding  to  those  of  Figure  8,  were: 

U  =  L2  =  0.333  ft 
L3  -  4.333  ft 
L  =  5.000  ft 
di  =  1.125  in 
d2  =  2.625  in 
d3  =  4.125  in 
d4  =  5.375  in. 


25 


Using  a  cast  specific  gravity  for  HBX-1  of  1.71  (Dobratz,  1981.  p.  19.53)  gives 
the  charge  weight  of  60.1  lbm.  Referring  to  Figure  17,  symmetric  boundary 
conditions  were  applied  on  the  yz  and  zx  planes.  Non-reflective  boundary 
conditions  were  applied  to  the  remaining  four  planes. 


WATER 


CHARGE 


Figure  17.  FEA  model  used  to  examine  the  three  dimensional  aspects  of 
the  tapered  charge. 

The  charge  and  water  model  were  analyzed  using  FEA  program 
VEC/DYNA3D  (Hallquist  and  Stillman.  1990).  This  is  the  stand-alone 
version  of  the  coupled  program  used  in  Chapter  III.  The  explosive  was 
modeled  as  a  high  explosive  burn  material  with  the  Jones-Wilks-Lee  (JWL) 
equation  of  state,  the  water  as  a  null  material  using  the  Gruneisen  equation 
of  state. 
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For  the  charge,  the  high  explosive  burn  material  model  used  by 
VEC/DYNA3D  requires  entry  of  the  detonation  velocity  D,  the  Chapman- 
Jouget  pressure  PCJ,  and  the  density  p  (Hallquist  and  Stillman,  1990,  p.  40). 
The  values  D  =  0.731  cm/us,  P(„  =  0.2204  Mbar,  and  p  =  1.712  gm/cm3  for 
HBX-1  were  taken  from  Dobratz  (1981,  p.  19.53).  The  JWL  equation  of  state 
was  used  to  describe  the  pressure-volume-energy  behavior  of  the  detonation 
products  (Dobratz,  1981,  p.  8.21): 

P  =  A(l-^)e-K.v+B(l-^)e-^  +  f 

Where  A,  B,  and  C  =  linear  coefficients  in  Mbar 

Ri ,  R2,  and  ©  =  nonlinear  coefficients 

v  volume  of  detonation  products 

v°      volume  of  undetonated  high  explosive 
P  =  pressure  in  Mbar 

E  =  detonation  energy  per  unit  volume  in . 

cm'3 

The  parameters  E,  A,  B,  Rr  R.„  and  co  listed  above,  empirically  derived  from 
cylinder-test  data,  are  required  entries  for  use  of  the  JWL  equation  of  state 
with  VEC/DYNA3D  (Hallquist  and  Stillman,  1990,  p.  89).  For  this  study 
these  values  were  taken  from  cylinder-test  data  for  H-6,  an  explosive  of 
similar  composition  to  HBX-1,  due  to  non- availability  of  HBX-1  data.  The 
values  used  were  (Dobratz,  1981,  p.  8.22): 

E  =  Q  103Mbarxjmii 

cm3 
A  =  7.5807  Mbar 
B  =  0.08513  Mbar 
Ri  =4.9 
R2  =  1.1 
co  =  0.20. 
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For  the  water,  the  VEC/DYNA3D  null  material  model  (Hallquist  and 
Stillman,  1990,  p.  41)  requires  entry  of  the  density  and  an  optional  pressure 
cutoff.  The  value  p  =  1.000  gm/crrr1  was  used  for  density,  and  a  pressure 
cutoff  value  of  6.89  x  10"9  Mbar  (0.1  psi)  were  used  since  water  is  unable  to 
sustain  tension.  The  Gruneisen  equation  of  state  with  cubic  shock  velocity 
-particle  velocity  (u -u  )  defines  pressure  p  in  Kbar  for  compressed  materials 
as  (Hallquist  and  Stillman,  1990,  p.  91): 


PoC2u 


1  +  11-7  V-lli2 


l-(S1-lJu-S2^-S3r^ 

Ui+1J 


+   yo  +  afi  E 


P 
where  li  =  - —  1 

^     Po 

i        •  gni 

p  =  density  in  — - 
p  cm3 

po  =  standard  density  in  — - 

cm3 

and  the  required  parameters  are 

C  =  the  intercept  of  the  us-up  curve  in  — 

Si,  S2,  S3=  coefficients  of  the  slope  of  the  us-up  curve 

Yo  =  the  Gruneisen  gamma 

a  =  first  order  volume  correction  to  yo . 

xr     ■    -         i                           +      1          •    Mbar  x  cm3 
E  =  internal  energy  per  unit  volume  in 

The  parameters  C,  Sr  S.„  and  S3,  for  the  us-up  curve  were  obtained  from 
Steinberg  (1987): 

us  =  C  +  SiUp+S2^up  +S3[^-J  Up 

=  0. 148  +  2.56up  -  1.986^-Up  +  0.2268^)  up 
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Appendix  B  outlines  the  determination  of  y,,  =  0.4934  and  a  =  1.3937,  based 
upon  the  seventh  order  polynomial  approximation  for  the  Gruneisen  gamma 
as  a  function  of  specific  volume  developed  by  Gurtman,  Kirsch.  and  Hastings 
(1971). 

B.     RESULTS 

The  large  number  of  elements  led  to  a  computationally  intensive  finite 
element  analysis  of  the  tapered  charge  underwater  explosion.  Running  on  a 
UNIX  engineering  workstation  required  approximately  four  days  to  perform 
the  computations  through  2  \xs  after  detonation.  Figure  18  shows 
representative  time  histories  of  two  elements  located  at  a  standoff  of  35"  from 
the  charge  nose.  One  element  is  on  the  charge  axis  off  the  nose,  the  other  90 
degrees  off  this  axis.  Due  to  the  large  variation  in  magnitude  of  pressures  at 
the  two  locations,  as  will  be  shown  in  Section  IV.B.2.  the  pressures  have  been 
normalized  to  better  compare  the  general  shapes  of  the  curves. 

The  element  located  on  the  charge  axis  maintained  a  greater  portion  of 
its  maximum  peak  pressure  over  a  longer  period  of  time  than  did  the  element 
located  90  degrees  off  the  axis.  This  observation  further  supports  the  use  as 
well  as  the  orientation  of  the  tapered  charge  for  simulation  of  the  underwater 
nuclear  explosion. 

Both  of  the  curves  of  Figure  18  rose  less  rapidly  than  expected  and 
displayed  oscillatory  behavior.  Gordon  and  Davidson  (1983)  experienced 
similar  results  when  analyzing  a  tapered  pentolite  charge  using 
finite -difference  techniques  in  two  dimensions.  They  attributed  the  long  rise 
time  to  three  possible  causes:  the  artificial  viscosity  coefficient  built  into 
their  model,  the  mesh  size,  and  the  fact  that  their  model  was  approximated 


29 


by  a  pointed  nose.  Although  the  model  used  in  this  study  did  not  use  a 
pointed  nose  approximation,  the  other  two  possible  causes  apply  to  this 
study.  Fox  (1992.  pp.  9.10)  cited  mesh  reflection  effects  as  a  contributing 
factor  to  instabilities  when  using  the  FEA  code  of  this  study. 
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Figure  18.  Normalized  Pressure-time  History  of  two  elements  located  at 
R  -  45  inches:  on  the  charge  axis  and  90  degrees  off  the  charge  axis. 


Although  desirable,  a  finer  mesh  was  beyond  the  capabilities  of  the 
machine  used  in  this  analysis.  For  this  study,  the  default  artificial  viscosity 
coefficient  and  time  integration  steps  (Hallquist  and  Stillman,  1990,  p.  25) 
were  used.  Future  sensitivity  studies  may  identify  values  for  these 
parameters.  The  mesh  reflection  effect  appears  to  be  a  strong  factor  in  the 
observed  pressure  oscillations.  The  key  contributor  to  this  problem,  uneven 
size  of  adjacent  elements,  was  aggravated  in  this  study  by  the  necessity  to 
conform  the  shape  of  eight-node  elements  to  the  curvatures  of  the  tapered 
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charge.  The  end  result  was  fine  charge  elements  near  the  nose  next  to 
relatively  coarse  water  elements  and  a  reversal  of  this  effect  at  the  tail  of  the 
charge. 

Oscillations  in  pressure  and  the  limited  time  of  the  pressure  histories 
generated  precluded  determination  of  the  pressure  plateau  duration.  The 
remainder  of  this  section  concentrates  on  the  shock  wave  as  it  emanates  from 
the  charge  at  early  times  and  the  directional  nature  of  the  peak  pressures  in 
the  media  surrounding  the  tapered  charge. 

1.      Early  Time  Shape  of  the  Shock  Front 

Figures  19  through  25  show  a  view  perpendicular  to  the  yz-plane 
of  the  water  elements.  The  model  has  been  reflected  about  the  zx-plane,  with 
charge  elements  removed.  Figures  19  applies  to  time  prior  to  detonation. 
Shading  indicates  the  elements  used  for  the  plots  of  Section  IV.B.2. 


Figure  19.   Water  elements  prior  to  detonation. 
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Figures  20  through  25  comprise  a  set  of  pressure  contour  plots  at 
early  times  through  0.58  ms.  There  are  five  contour  lines  in  each  plot.  The 
contours  range  from  1.000  psi  for  the  sparsest  dotted  line  to  5,000  psi  for  the 
solid  line.  This  relatively  low  pressure  range  at  such  close  proximity  to  the 
charge  was  selected  to  clearly  define  the  location  of  the  evolving  shock. 


Figure  20.   Pressure  contour  plot  of  water  elements  at  t  =  0.08  ms  after 
detonation. 


The  effect  of  mesh  reflection,  discussed  in  the  previous  section,  can 
be  seen  in  Figure  20  along  the  charge  axis  to  the  front  of  the  charge  nose  as 
the  shock  wave  transmits  through  the  media  at  different  speeds  in  this 
region  of  very  poor  match  of  element  sizes.  The  pressure  contours  bend 
inward  toward  the  charge  nose  in  the  small  element  region.  By  0.18  ms, 
Figure  21,  the  mesh  reflection  effect  in  front  of  the  charge  is  negligible. 
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Figure  21.   Pressure  contour  plot  of  water  elements  at  t  =  0.18  ms  after 
detonation. 


Figure  22.   Pressure  contour  plot  of  water  elements  at  t  =  0.28  ms  after 
detonation. 
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Figure  23.   Pressure  contour  plot  of  water  elements  at  t  -  0.38  ms  after 
detonation. 
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Figure  24.  Pressure  contour  plot  of  water  elements  at  t  =  0.48  ms  after 
detonation. 
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Figrure  25.   Pressure  contour  plot  of  water  elements  at  t  =  0.58  ms  after 
detonation. 


The  shock  front  begins  as  a  tear  shape  at  early  times  emanating 
from  the  charge  burn  region  and  becoming  nearly  spherical  on  the  charge 
axis  off  the  nose.  The  above  pressure  contours  show  the  nose  portion  of  the 
evolving  shock  front  to  be  nearly  spherical  in  shape  out  to  approximately  40 
degrees  off  the  charge  axis  with  a  radius  centered  at  the  location  of  the 
undetonated  charge  nose.  Once  the  burn  region  reaches  the  tail  of  the 
charge,  the  aft  portion  of  the  shock  front  expands  to  nearly  spherical  with  a 
smaller  radius  of  curvature  than  at  the  forward  end.  The  above  figures 
clearly  illustrate  the  shock  rapidly  travelling  outward  ahead  of  the 
expanding  detonation  products. 
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2.      Directionality  of  Peak  Pressure 

This  subsection  examines  variations  in  the  peak  pressure 
developed  in  regions  surrounding  the  tapered  charge  based  upon  the  FEA 
conducted.  Figure  19  marks  the  locations  of  the  elements  whose 
pressure-time  histories  were  gathered  and  processed  to  form  the  plots. 

For  the  plots  of  this  section,  relative  peak  pressure  Pre]  is  defined 
as: 


PrPi  = 


Po 


Where 


Pmax  =  maximum  pressure  computed  for  the  element 
P0  =  maximum  pressure  observed  for  the  element  at  the  same 
standoff  and  nearest  the  charge  nose  axis. 


Figures  accompanying  the  plots  serve  to  further  clarify  the  determination  of 
Pre!  as  well  as  explain  the  geometry  corresponding  to  each  plot. 

Figure  26  corresponds  to  Figure  27,  a  plot  of  relative  peak  pressure 
as  a  function  of  9  degrees  off  the  charge  nose  axis  from  9  =  0  to  90  degrees,  at 
constant  standoff,  for  three  different  standoffs. 


<-- 


Constant  R 


Target 


Figure  26.  Constant  standoff,  6  degrees  off  charge  axis. 
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Figure  27.    Relative  peak  pressure  as  a  function  of  the  angle  off  the 
charge  axis  at  constant  standoff. 


For  the  three  standoff  distances  of  Figure  27,  there  is  a  decrease  in 
pressure  from  the  on-axis  pressure  as  the  off-axis  angle  increases  from  0  to 
20  degrees.  This  decrease,  ranging  from  3  per  cent  for  R  =  43  inches  to  5  per 
cent  for  R  =  15  inches,  may  be  attributed  to  mesh  geometry  and  mesh 
reflections.  The  effect  of  mesh  geometry  was  caused  by  the  centroid  of  each 
element  not  being  at  the  exact  standoff.  This  effect  was  minimized  by 
applying  a  first  order  correction  to  each  measured  pressure: 

Where  Pr  =  pressure  used  for  plot  at  the  nominal  standoff  R 

PR'  =  max  pressure  for  the  given  element  at  standoff  R 

R  =  Jx2  +y2+z2 

x,  y,  z  =  coordinates  of  the  element  centroid. 
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After  approximately  20  degrees,  the  curves  of  Figure  27  rise 
steadily  with  increasing  angle  off  the  charge  axis.  This  was  to  be  expected 
due  to  the  tradeoff  between  peak  pressure  duration  and  magnitude  of  the 
tapered  charge.  As  0  increases,  the  target  is  placed  nearer  to  the  bulk  of  the 
mass  of  the  charge.  In  a  test  situation  gages  must  often  be  located  off  the 
charge  axis  to  record  the  free-field  pressure  experienced  by  the  target 
without  being  overly  influenced  by  proximity  to  the  target.  The  test  designer 
often  attempts  to  locate  the  gage  as  far  off  axis  as  possible  while  minimizing 
the  deviation  of  the  gage  measurement  from  the  on-axis  values.  Figure  27 
indicate  that,  for  the  charge  and  standoffs  of  this  study,  gages  could  be 
located  up  to  35  degrees  off  the  charge  axis  for  an  error  in  maximum  peak 
pressure  measurement  of  less  than  5  per  cent.  To  keep  the  error  in 
measurement  of  plateau  duration  to  an  acceptable  level,  the  allowable 
off-axis  angle  may  lie  well  below  the  level  based  upon  peak  pressure  alone. 

Past  90  degrees  off  the  front  charge  axis,  constant  standoff  was 
replaced  by  constant  distance  y  off  the  charge  axis  to  examine  the  relative 
pressures  experienced  by  elements  to  the  side  of  the  charge.  Figure  28 
illustrates  the  geometry  involved.  Po  for  this  case  is  taken  at  a  standoff  R  =  y 
on  the  axis  off  the  charge  nose  as  shown  in  the  figure.  Figure  29  shows 
relative  peak  pressure  as  a  function  of  £,  /  L  off  the  charge  axis  for  \  I  L  =  0  to 
1  from  nose  to  tail.  As  can  be  seen  from  the  figure,  relative  peak  pressure 
shows  marked  departure  for  different  values  of  y.  While  the  highest  relative 
peak  pressure  occurred  at  £,  /  L  =  0.8,  the  magnitude  fell  from  8.4  at  y  =  19 
inches  to  6.2  at  y  =  43  inches.  This  was  likely  a  close-in  phenomenon.    It  is 
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Figure  28.     Constant  distance  off  the  charge  axis,  to  the  side  of  the 
charge. 
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Figure  29.    Relative  peak  pressure  as  a  function  of  the  fraction  of  the 
length  aft  of  the  nose  at  constant  distance  off  the  charge  axis  to  the  side. 


expected  that,  at  distances  further  removed  from  the  charge  axis  than  those 
of  this  study,  the  maximum  relative  peak  pressure  to  the  side  of  the  charge 
will  continue  to  fall  then  reach  a  steady  value  significantly  lower  than 
calculated  here.     Of  interest  is  the  fact  that  the  maximum  relative  peak 
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pressure  observed  for  elements  located  to  the  side  of  the  charge  was  found  at 
£  /  L  =  0.8,  significantly  aft  of  the  center  of  gravity  location  at  £  /  L  =  0.6. 
This  may  be  due  to  a  build-up  in  pressure  away  from  the  center  of  gravity 
location  in  the  direction  of  the  tapered  charge  high  explosive  burn. 

Completing  the  examination  of  pressures  encircling  the  charge 
from  fore  to  aft,  Figure  30  shows  the  geometry  corresponding  to  Figure  31 
plotting  the  variation  in  relative  peak  pressure  at  a  constant  standoff  R 
measured  from  the  center  of  the  tail.  The  angle  0  in  this  case  is  measured 
from  90  degrees  off  the  rear  charge  axis  toward  the  rear  charge  axis  in  order 
to  continue  proceeding  in  a  counterclockwise  direction.  As  in  the  two  other 
cases,  Po  is  measured  off  the  charge  nose,  in  this  case  at  a  standoff  equal  to 
the  constant  R.  As  in  the  constant  R  case  off  the  charge  nose,  a  correction 
was  applied  to  the  calculated  values  to  account  for  variation  in  standoff  of  the 
element  centroids. 
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Figure  30.  Constant  R  off  the  rear  charge  axis. 


Figure  31  shows  a  continuation  of  the  decrease  in  peak  pressure 
for  the  standoffs  studied  until  0  =  60  degrees,  or  30  degrees  off  the  charge  tail 
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Figure  31.    Relative  peak  pressure  as  a  function  of  degrees  aft  of  the 
charge. 

axis  followed  by  a  slight  increase  to  the  charge  axis.  Here  again,  the 
unevenness  of  the  curves  near  the  axis  may  be  attributed  in  part  to  mesh 
reflections.  As  in  the  side  case,  the  relative  pressure  curves,  though 
exhibiting  similar  shape,  varried  significantly  for  different  standoffs,  with 
shorter  standoff  corresponding  to  higher  relative  peak  pressure.  Though  of 
minor  interest  in  nuclear  simmulations,  the  side  and  rear  relative  pressure 
plots  are  of  more  interest  to  weapons  and  industrial  high  explosive  designers. 
The  higher  pressures  at  the  tail  end  of  the  charge  agree  with  the  use  of  this 
configuration  for  demolition  and  other  applications. 

Of  more  interest  in  simulations  is  the  effect  of  distance  off  axis  at  a 
constant  standoff  from  the  charge  nose  as  illustrated  in  Figure  32.  This 
configuration,  somewhat  easier  to  set  up  for  experimental  verification, 
yielded  the  plot  of  Figure  33  which  provides  information  of  similar  utility  as 
the  constant  R  plot  of  Figure  27. 
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Figure  32.  Constant  distance  from  nose  along  charge  axis. 
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Figure  33.    Relative  peak  pressure  at  constant  distance  along  charge 
axis  as  a  function  of  degrees  off  the  charge  axis. 


Figure  33  indicates  that,  for  the  studied  charge  design  and  stand 
off  distances,  the  off  axis  distance  may  be  varied  up  to  10  degrees  before 
exceeding  a  5  per  cent  error  in  peak  pressure  readings.  Since  the  allowable 
angle  to  stay  within  a  given  error  tolerence  decreased  with  range,  a  smaller 
angle  would  be  allowed  at  the  longer  ranges  used  in  simulations  of  the 
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underwater  nuclear  explosion.  Having  concluded  the  FEA  of  the  underwater 
tapered  charge  explosion,  the  next  chapter  outlines  computer  optimization  of 
a  simple  method  to  determine  the  pressure-time  history  of  a  tapered  charge. 
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V.  TAPERED  CHARGE  DESIGN  OPTIMIZATION 

This  chapter  outlines  the  coupling  of  a  public  domain  computer 
optimization  package  ADS  (Vanderplaats.  1984)  and  a  simple  tapered  charge 
pressure-time  generating  program  TAPER  (Costanzo,  1991)  to  optimize  the 
tapered  charge  design  process.  The  coupling  program,  listed  in  its  entirety  in 
Appendix  C,  can  be  used  with  existing  and  future  tapered  charge 
pressure-time  history  generating  codes  to  enhance  the  design  of  the  tapered 
charges  used  to  simulate  underwater  nuclear  explosions. 

In  addition  to  the  program  in  Appendix  C.  ADS  (Automated  Design 
Synthesis)  and  a  pressure-time  generating  subroutine  are  required  to 
perform  tapered  charge  design  optimization. 

A.     SIMPLE  PRESSURE-TIME  HISTORY  ALGORITHM 

Simpler,  less  computationally  intensive,  computer  codes  than  that  used 
for  the  finite  element  analysis  of  the  previous  chapter  exist  to  predict  the 
pressure  profile  of  an  underwater  tapered  charge  explosion.  These  simpler 
codes  are  usually  based  upon  an  empirically  based  superposition  principle. 

Because  the  empirically  based  exponential  approximation  of  Chapter  II 
works  well  for  spherical  /  near  spherical  charges,  one  method  of  deriving  the 
pressure-time  history  of  an  underwater  tapered  charge  explosion  is  to 
partition  the  tapered  charge  into  subsegments  as  shown  in  Figure  34.  Each 
subsegment  is  considered  to  be  a  separate  charge  generating  its  own 
pressure-time  history.  (Costanzo.  1991) 
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Figure  34.    Tapered  charge  discretized  into  subsegments  for  calculation  of 
pressure-time  history  using  summation  method. 


Shocks  generated  by  subsegments  detonated  after  the  nose  subsegment 
travel  in  shocked  media  at  faster  speeds  than  preceding  shocks,  overtaking 
them.  There  is  thus  a  summation  or  stacking  effect  of  individual  "wavelets" 
to  determine  the  overall  pressure-time  history  of  the  tapered  charge 
underwater  explosion  (Costanzo.  1991): 


£  A- 


1=1 


ri 


S£ 


Where 


P  =  pressure  as  a  function  of  time 

N  =  number  of  waves  stacked 

A,  =  empirical  amplification  factor,  a  function  of  0, 

9,  =  slope  angle  of  subsegment 

r,  =  R  +  £,,  =  subsegment  standoff 

R  =  standoff  from  charge  nose  to  target 

£,  =  distance  from  nose  to  subsegment  midpoint 

B  =  empirical  decay  constant 

t  =  time 

Dj  =  charge  diameter  at  subsegment  midpoint 

8£  =  subsegment  length. 

Figure  35  illustrates  the  superposition  scheme  for  a  simple  tapered  charge 
discretized  into  large  segments. 
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Figure  35.   Summation  of  wavelets  generated  by  charge  subsegments  to 
obtain  resultant  pressure-time  history  for  a  tapered  charge. 

One  algorithm,  using  a  superposition  scheme  and  empirical  data  to 
determine  the  pressure-time  history  of  a  tapered  charge  detonated 
underwater  is  the  PASCAL  program  TAPER  written  by  Fred  Costanzo 
(1991).  This  program  was  converted  to  a  FORTRAN  program  then  to  the 
FORTRAN  subroutine  used  for  tapered  charge  design  optimization.  Besides 
a  basic  overhaul  of  the  input  /  output,  a  simple  provision  for  surface  cutout 
was  added  to  the  original  algorithm. 

B.     OPTIMIZATION  PROBLEM 

The  objective  of  the  tapered  charge  design  optimization  was  to 
determine  charge  geometry,  as  depicted  in  Figure  8,  and  standoff  required  to 
develop  a  pressure-time  history  most  nearly  matching  a  desired  pressure- 
time  history.  The  average  square  root  of  the  sum  of  the  squares  of  the 
differences   between   the   computed   and   desired  pressures    at   each   time 
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increment  was  selected  to  quantify  the  difference  between  the  optimal  design 
and  desired  pressure  profiles.  Constraints  restricted  the  charge  to  increasing 
diameters  and  placed  upper  and  lower  limits  on  charge  weight,  charge 
dimensions,  and  standoff.  Only  designs  within  these  limits  were  allowed. 
For  a  tapered  charge  of  three  segments,  the  optimization  problem  became: 


minimize:       ,  2L, 77 >  i  =  1,  N 


subject  to:      1  in  <  L  <  20  ft,  i  =  1  to  3 

0.75  in  <  d  <  10  in,  i  =  1  to  4 

5  ft  <  R  <  30  ft, 
25  lbm  <  W  <  125  lbm. 
d,  <  d,  <  dr  <  d, 

where  N  =  number  of  computed  pressure-time  history  points 

Pj  =  desired  pressure  at  a  particular  time 
p  =  computed  pressure  at  a  particular  time 
L  =  charge  segment  length 
di  =  charge  joint  diameter 
R  =  standoff 
W  =  charge  weight. 

The  main  FORTRAN  program  DTAPOPT  was  written  to  accomplish  the 
tapered  charge  design  optimization  using,  as  described  previously,  the  ADS 
optimization  package  and  the  subroutine  based  upon  the  TAPER  program. 

C.     OPTIMIZATION  PROGRAM 

In  addition  to  DTAPOPT,  three  subroutines  for  use  in  conjunction  with 
the  main  program  were  also  written.  The  first.  DCOMPAR.  computes  the 
square  root  of  the  average  sum  of  the  squares  of  the  pressure  differences. 
The     second,     DPTGEN,     interpolates     to     find     the     desired     pressure 
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corresponding  to  the  design  pressure  at  a  given  time.     This  subroutine 
enables  the  comparison  of  the  desired  and  design  pressures  at  the  same 
times.    The  third  subroutine,  DTAPWT,  computes  the  weight  of  the  tapered 
charge  for  evaluation  of  the  weight  constraint. 
1.      Required  Program  Input 

DTAPOPT  requires  input  from  two  files  and  the  keyboard.  Figure 
36  shows  a  sample  of  a  tapin.dat  file  used  to  input  initial  charge  geometry, 
standoff,  nominal  length  of  pressure-time  history  to  be  computed,  and  surface 
cutoff  time.  Decimal  alignment  and  horizontal  placement  on  the  lines  is 
optional  with  one  number  per  line  only.  Number  translations  have  been 
written  into  the  figure. 


3 

number  of  charge  segments 

0.5 

length  of  first  segment  (ft) 

1.0 

length  of  second  segment  (ft) 

5.0 

length  of  third  segment  (ft) 

1.0 

first  joint,  or  nose,  diameter  (in) 

2.0 

second  joint  diameter  (in) 

3.0 

third  joint  diameter  (in) 

5.0 

fourth  joint,  tail,  diameter  (in) 

15.1 

standoff  (ft) 

1.12 

nominal  length  of  computed  time  history  (ms) 

1.095 

surface  cutoff  time  (ms) 

Figure  36.    Sample  tapin.dat  file.   Italics,  not  part  of  the  file,  indicate 
the  meaning  of  each  number. 

Figure  37  shows  a  sample  of  the  second  input  file,  profin.dat,  used 
to  input  the  desired  pressure  profile  for  program  DTAPOPT.  Up  to  1001  data 
points  may  be  entered  without  program  modification.  The  first  time  entered 
must  be  at  time  zero,  and  the  last  time  must  be  greater  than  the  nominal 
length  of  computed  time  history  entered  in  tapin.dat.     Decimal  alignment 
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and  horizontal  placement  on  the  lines  is  optional.  Only  one  number  can  be 
used  on  the  first  line,  two  on  the  rest.  Number  translations  added  to  the 
figure  are  shown  in  italics. 


tides,  total  number  of  desired,  pressure-time  data  points  minus  1 


0.0             0.0 

tdes(O)  in  ms--musi  be  zero. 

pdes(0)  in  psi 

0.04     1950.0 

tdes(l) 

pdes(J) 

0.119  1755.0 

tdes(2) 

pdes(2) 

0.278  1560.0 

tdes(S) 

pdes(3) 

0.417  1365.0 

tdes(4) 

pdes(4) 

0.635  1170.0 

ldes(5) 

pdes(5) 

0.834     975.0 

tdes(6) 

pdes(6) 

1.072     780.0 

tdes(7) 

pdes(7) 

1.120         0.0 

tdes(H) 

pdes(8) 

3.0             0.0 

tdes(9) 

pdes(9) 

tdes(ndes)  must  In' greater  than  nominal  length  ofp-t  history  entered  in  tapin.dat 


Figure  37.    Portion  of  a  sample  profin.dat  file.    Italics,  not  part  of  the 
file,  have  been  added  to  indicate  what  each  number  represents. 


Three  integers  comprise  required  keyboard  input.  These  numbers 
control  the  optimization  method  used  by  ADS.  Vanderplaats  (1985)  provides 
more  detailed  instruction  on  method  selection.  Vanderplaats  (1984)  provides 
the  theory  behind  the  methods.  The  three  numbers  consist  of  any  one 
number  from  each  of  three  groups.  The  first  number  may  be  any  of  the 
following  to  determining  the  optimization  strategy: 


First 
Number  Optimization  Strategy 

0  Go  directly  to  the  optimizer 

6  Sequential  Linear  Programming 

7  Method  of  Centers  (Design  must  be  feasible) 

8  Sequential  Quadratic  Programming 

9  Sequential  Convex  Programming. 
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The   second   required   input   number   may   be   either   of  the   following   to 
determine  the  optimizer: 

Second 
Number  Optimizer 

4  Method  of  Feasible  Directions 

5  Modified  Method  of  Feasible  Directions. 

The  last  input  number,  one  of  the  following,  selects  the  one-dimensional 
search  option  to  be  used: 

Third 
Number  One-Dimensional  Search 

5  Golden  Section  Method 

6  Golden  Section  Method  Plus 

Polynomial  Interpolation 

7  Bounded  Polynomial  Interpolation 

8  Unbounded  Polynomial  Interpolation. 

2.      Program  Output 

Output  from  DTAPOPT  includes  one  screen  summarizing  the 
optimization  and  an  output  file  containing  the  design  and  interpolated 
desired  pressure-time  histories. 

Figure  38  shows  a  sample  screen  from  a  run  of  DTAPOPT  on  a 
personal  computer.  The  figure  shows  the  execution  commands  followed  by 
prompts  for  the  three  inputs,  in  this  case  the  user  selected  the  combination 
8-5-7  for  the  strategy,  optimizer,  and  search.  The  output  summary  then  lists 
the  initial  and  final  charge  designs  as  well  as  the  average  square  root  of  the 
pressure  differences  for  each  case  and  the  number  of  calls  to  the  pressure- 
time  generating  routine. 
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C : \P&C\FOR\OPT : dtapopt 

ENTER  ISTRAT.IOPTJONED:   8  5  7 

INITIAL  DESIGN 

LENGTHS  =        1.000     1.000 

5.000 

FT 

DIAMETERS  -       1.00      2.00 

4 . 00      5 . 00 

IN 

CHARGE  WEIGHT  =  66.0  LB 

RANGE  =         20.0  FT 

SQRT  OF  AVG  OF  (Pdesign-Pdesired)A2 

=  424.3  PSI 

FINAL  DESIGN 

LENGTHS  =        0.083     0.083 

4.825 

FT 

DIAMETERS  =      2.52      3.75 

4.28      4.52 

IN 

CHARGE  WEIGHT  =  55.7  LB 
RANGE  =  17.4  FT 
SORT  OF  AVG  OF  (Pdesign-Pdesired)A2  =  240.3  PSI 


CALLS  TO  P-T  HISTORY  GENERATOR 
C:\P&C\FOR\OPT: 


158 


J 


Figure  38.  Sample  screen  output  from  DTAPOPT. 


The  output  screen  of  Figure  38  represents  an  early  step  in  the 
design  process.  The  next  would  be  to  use  this  final  design  to  input  a  more 
refined  initial  design,  then  run  the  optimization  program  again.  For  test 
optimizations,  the  square  root  of  the  average  pressure  difference  squared  was 
well  below  100  for  the  "best"  final  design. 

Figure  39  shows  a  portion  of  a  sample  DTAPOPT  output  file.  The 
first  two  columns  contain  the  calculated  time  and  pressure,  the  third  column 
contains  the  interpolated  values  of  the  desired  pressure  profile  input.  The 
data  are  in  free  format  to  retain  maximum  precision,  sacrificing  the 
readability  of  formatting.  This  output  file  may  be  readily  used  to  create  plots 
comparing  the  initial  and  final  designs  as  was  done  in  the  next  section  of  this 
chapter. 
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0 . OOOOOOOOOOOOOOOOOE -010. OOOOOOOOOOOOOOOOOE - 01  0 . OOOOOOOOOOOOOOOOOE - 01 
2 . 81022269002981981E-02  1639 . 32749805998742000  1405 . 11134501490983000 
5 . 62044538005963962E-02  1675 . 89437624031484000  1987 . 35262142392480000 


0 . 47773785730506935  1628 . 10858635031309000  1658 . 35094063994575000 
0 . 50584008420536752  1614 . 83610147863396000  1636 . 41749525434739000 
0 . 53394231110566581  1602 . 96393658242482000  1614 . 48404986874857000 


1.06788462221133162   900.78633515776994000   899.17511224684437800 
1 . 09598684911162980  0 . OOOOOOOOOOOOOOOOOE -01   235 . 12876911529235700 


Figure  39.    Sample  excerpted  from  an  output  data  file  generated  by 
DTAPOPT. 


D.     OPTIMIZATION  RESULTS 

The  design  space  using  the  subroutine  adapted  from  program  TAPER 
proved  to  be  fraught  with  local  minima  attributable  discontinuities  resulting 
from  integer  changes  in  the  number  of  waves  stacked.  Several  runs  of 
DTAPOPT  using  different  optimization  methods  for  a  particular  initial 
design  resulted  in  an  improved  design  unless  the  initial  design  was  optimal. 
The  improved  design  was  then  used  as  the  initial  design  for  further 
optimization.  This  process  was  repeated  from  four  to  seven  times  until 
further  optimizations  failed  to  improve  the  design  an  appreciable  amount. 
Each  run  of  DTAPOPT  took  approximately  two  minutes  on  a  personal 
computer  of  modest,  386SX,  capacity.  Total  time  to  perform  a  tapered  charge 
design  optimization  was  from  one  to  two  hours. 

Of  the  desired  pressure  profiles  and  initial  designs  tested,  input 
parameter  combinations  0-5-7,  8-5-7,  and  9-5-7,  usually  produced  the  most 
improved  design  with  the  fewest  calls  to  the  pressure-time  history  generating 
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subroutine.  These  input  numbers  translate,  as  shown  in  Section  V.C.I  to 
direct  optimization,  sequential  quadratic  programming,  or  sequential  convex 
programming  strategy  combined  with  the  modified  method  of  feasible 
directions  optimizer  and  a  bounded  polynomial  interpolation 
one-dimensional  search.  For  the  best  designs  found,  as  mentioned  in  Section 
V.C.2,  the  square  root  of  the  average  difference  between  desired  and  design 
pressures  was  well  below  100  psi.  Figure  40  shows  a  comparison  of  the 
pressure  profile  generated  by  an  optimized  design  from  program  DTAPOPT 
and  the  corresponding  desired  pressure  profile. 
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Figure  40.   Pressure  profile  resulting  from  optimization  compared  with 
the  corresponding  desired  pressure  profile. 


53 


VI.  CONCLUSIONS  AND  RECOMMENDATIONS 

The  boundary  element  /  finite  element  analysis  of  the  side  on  attack  of  a 
simple  cylinder  by  nuclear  and  conventional  type  pressure  profiles  resulted 
in  substantially  different  structural  responses.  It  is  therefore  necessary  to 
model  the  duration  of  the  pressure  plateau  to  adequately  simulate  structural 
response  to  an  underwater  nuclear  explosion. 

The  finite  element  analysis  of  an  underwater  tapered  charge  explosion 
provided  insight  to  the  early  time  propagation  of  the  shock  wave  and  the 
directional  variations  in  peak  pressures  developed  in  the  surrounding  media. 
There  is  need  for  future  research,  including  sensitivity  studies  of  the 
artificial  viscosity  coefficient  and  the  time  integration  step,  to  obtain  less 
oscillatory  results  to  provide  tapered  charge  pressure  plateau  duration 
information.  Additionally,  the  computationally  intensive  method  used  lends 
itself  more  readily  to  supercomputer  use  where  the  mesh  size  may  be 
extended  far  enough  away  from  the  charge  form  comparison  with  test  data. 

Coupling  of  an  optimization  routine  with  a  tapered  charge  pressure 
profile  generating  routine  provides  a  tool  which  can  be  used  to  more 
efficiently  design  tapered  charges  used  to  simulate  underwater  nuclear 
explosions.  The  program  written  for  this  study  may  be  used  with  existing 
and  future  pressure-time  history  codes. 


54 


APPENDIX  A:  TAPERED  CHARGE  FEM  INPUT  FILE 

Following  is  a  complete  listing  of  the  input  file  used  to  perform  FEM 
analysis  on  the  tapered  charge  and  water  model  of  Chapter  IV: 


GOG02 

c  Set  up  for  DYNA3D.  integrate  to  time  2000,  TAURUS  data  dump  interval 

10. 

c   high  speed  printer  dump  interval  29999. 

dn3d  vec  term  2000  plti  10.0  prti  29999.0 

c  Scale  x,  y.  z  axes  from  inches  to  cm. 

xsca  2.54  ysca  2.54  zsca  2.54 

c  Define  material  1,  type  8- -High  explosive  burn.  HBX-1: 

c   detonation  vel  0.731  cm/us.  Chapman- Jouget  pres  0.2204  Mbar . 

c   density  1.712  gm/cc. 

mat  1  type  8  d  .731  pcj  .2204  ro  1.712 

c  Equation  of  state  for  charge:  2--JWL: 

c   a  7.5807.  b  .08513.  rl  4.90.  r2  1.10.  omega  .20.  eO  .103  Mbar-cc/cc. 

eos  2  a  7.5807  b  .08513  rl  4.90  r2  1.10  omega  0.20  eO  0.103  endmat 

c  Define  material  2,  type  9- -Null  Material.  Water: 

c   cutout  pressure  .1  psi  (6.89e-9  Mbar).  density  (3  4  deg  C  1.000  gm/cc. 

mat  2  type  9  pc  6.89e-9  ro  1.000 

c  Equation  of  state  for  water:  4- -Gruneisen: 

c   sound  speed  .142  cm/us  (9  4  deg  C.  si  2.56.  s2  -1.986.  s3  .2268. 

c   Gruneisen  gamma  .4934,  first  order  vol  cor'n  to  gamma  1.3937. 

eos  4  sp  .142  si  2.56  s2  -1.986  s3  0.2268  gamma  .4934  sa  1.3937  endmat 

c  Define  two  symmetry  planes  by  defining  a  point  and  a  normal  vector 

c   for  each  plane.   Any  point  within  .001  cm  of  a  symmetry  plane  will  be 

c   included  in  that  plane's  definition. 

plane  2 
0  0  0   1  0  0  .001  svmm 
0  0  0   0  10  .001  symm 

c  Define  detonation  points  in  HBX-1. 

detp  1  point  0  0  0: 
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c  Define  rear  and  front  water  cylinders,  cone  and  plane  surfaces  for 
c   charge  nose  (part  1)  and  transition  water  (part  2). 


sd  1  cn2p  0  0  0  0  0  1 
sd  3  plan  0  0  0  0  0  1 
sd  4  plan  0  0-4   0  0  1 


0.5625 


0.0   1.3125 


■4.0 


c  Part  1:  small  charge  cone. 


start 
13  5; 
13  5: 
1  3; 
-10  1 
-10  1 
0.0    -4.0 

di    1   0   3;    10   3: 
sfi      -1    -3:    -1    -3 
sfi      -1    -3:    -1    -3 
sfi      -1    -3:    -1    -3 
dill      232 
dill      322 
mate    1 

end 


sd  1 
■1  -1;  sd  3 
■2  -2:  sd  4 


c  Part  2:   water  transition  from  small  charge  cone  to  square  grid. 


start 
13  5  7  9: 
13  5  7  9: 
1  3: 

-4-4044 
-4-4044 
0.0  -4.0 

di  12045:12045 
d  2  2  0  4  4  0 


sfi      -2  -4 

sfi      -1  -5 

sfi      -1  -5 

dill  352 

dill  532 
mate    2 
end 


2  -4 
1  -5 
1  -5 


;  sd  1 
-1  -1 
-2  -2 


sd  3 
sd  4 


c  Define  cone  and  plane  surfs  for  med  charge  and  transition  (pts  3  and 
4). 


sd  8   plan  0  0-8   0  0  1 

sd  9   cn2p  0  0   0   0  0  1    1.3125 

c  Part  3:  medium  charge  cone. 

start 
13  5; 
13  5; 
1  3; 
-10  1 
-10  1 
-4.0  -8.0 


■4.0    2.0625 


8.0 
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di    1   0   3:    1   0   3; 


sfi      -1    -3 
sfi      -1    -3 
sfi      -1    -3 
dill 
dill 
mate    1 
end 


-1  -3 
-1  -3 
-1  -3 

2  3  2 

3  2  2 


:  sd  9 

-1  -1:  sd  4 

-2  -2:  sd  8 


c  Part  4:  medium  water  transition. 

start   13    5    7   9: 
13    5    7    9; 
1    3; 

-4-4044 
-4-4044 
-4.0    -8.0 

di    12045:12045;    ; 
d   2    2   0     4  4   0 


sfi   -2 

-4  :  -2 

-4  : 

:  sd  9 

sfi   -1 

-5  :  -1 

-5  : 

-1  -1  : 

sd  4 

sfi   -1 

-5  :  -1 

-5  : 

-2  -2  : 

sd  8 

dill 

3  5  2 

dill 

5  3  2 

mate  2 

end 

c   Define    cones   and   plane   for   large    charge   and   transition    (pts    5   and    6) 

2.0625  -8.0        2.6875  -60.0 


sd  11  cn2p  0  0 

0 

0  0  1 

sd  13  plan  0  0 

-60 

0  0  1 

c  Part  5:  large  charge  cone. 


start 
13    5: 
13   5: 
1   27; 
-10   1 
-10   1 
-8.0    -60.0 
di   1   0   3:    10   3: 
sfi      -1    -3;    -1    -3 
sfi      -1    -3;    -1    -3 
sfi      -1    -3:    -1    -3 
dill      232 
dill      322 
mate    1 

end 


;  sd  11 

-1  -1;  sd  8 

-2  -2:  sd  13 


c  Part  6:  large  water  transition. 

start 
13    5   7   9: 
13   5   7   9; 
1    27: 

-4-4044 
-4-4044 
-8.0    -60.0 
di   12045:12045:    ; 
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d   2    2   0     4  4   0 

sfi      -2    -4    :    -2    -4    :     ;    sd   11 
sfi      -1    -5    ;     -1    -5    ;    -1    -1    ;    sd   8 
sfi      -1    -5    :    -1    -5    :    -2    -2    ;    sd   13 
dill      352 
dill      532 
mate    2 
end 

c  Define  cylinders  for  front  and  rear  water  cylinders  and  transitions. 

sd  cl  cyli  0  0  0   0  0  1    2.6875 

c  Part  7:  rear  water  cylinder. 

start 

13    5; 

13   5: 

1    24; 

-10   1 

-10   1 

-60.0    -106.0 

di    1   0   3:    1   0   3;     ; 

sfi      -1    -3:    -1    -3:    ;    sd  cl 

sfi      -1    -3:    -1    -3:    -1    -1:    sd   13 

dill      232 

dill      322 

nrb   0   0    2      0   0    2 

mate    2 
end 

c  Part  8:  rear  water  transition. 

start 

13   5   7   9: 

13    5   7   9: 

1    24: 

-4-4044 

-4-4044 

-60.0    -106.0 

di    12045:    12045;     ; 

d   2    2   0  4   4   0 

sf    2    2    1      4   4    2      sd   cl 

sfi    -1    -5;    -1    -5:    -1    -1;    sd   13 

d      111      352 

d      111      532 

nrb   0   0    2      0   0    2 

mate    2 
end 

c  Define  cones  and  plane  for  water  cone  and  transition  (parts  9  and  10) 

sd  5  cn2p  0  0   0    0  0  1    0.5625       0.0    1.3125        4.0 
sd  7  plan  0  0   4    0  0  1 

c  Part  9:  front  water  cone. 

start 
13  5; 
1  3  5; 
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1    3: 
-10   1 
-10   1 

0.0  4.0 
di    1   0   3: 
sfi      -1    ■ 
sfi      -1    - 
sfi      -1    ■ 
dill 
dill 
mate    2 
end 


10  3: 
3:  -1  -3 
3;  -1  -3 
3;  -1  -3 

2  3  2 

3  2  2 


sd  5 

1  -1:  sd  3 

2  -2;  sd  7 


c  Part  10:  front  water  cone  transition. 

start 
13    5   7   9: 
13   5   7   9: 
1   3: 

-4-4044 
-4-4044 
0.0   4.0 

di   12045:12045:    : 
d   2    2   0     4  4   0 


sfi      -2    -4 
sfi      -1    -5 
sfi      -1    -5 
dill 
dill 
mate    2 
end 


-2  • 

-1 

-1  • 
3  5  2 
5  3  2 


:  sd  5 

-1  -1:  sd  3 

-2  -2:  sd  7 


c  Define  cylinder  for  front  water,  parts  11  and  12 
sd  c2  cyli  0  0  0   0  0  1    1.3125 
c  Part  11:  front  water  cylinder. 


start 
13    5: 
13    5; 
1   22: 
-10   1 
-10    1 
4.0   46.0 
di    1   0   3;    1   0   3: 
sfi    -1    -3;    -1    -3 
sfi    -1    -3;    -1    -3 
dill      232 
dill      322 
nrb  0   0   2      0   0   2 
mate    2 

end 


sd  c2 
1  -1:  sd  7 


c  Part  12:  front  water  cylinder  transition 

start 
13  5  7  9; 
13  5  7  9; 
1  22: 
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-4-4044 
-4-4044 
4.0   46.0 

di   12045:12045:    ; 
d   2    2   0  4  4   0 
sf    2    2    1      4   4    2      sd   c2 
sfi    -1    -5:    -1    -5:    -1    -1:    sd   7 
d      111      352 
d      111      532 
nrb  0   0   2      0   0   2 
mate    2 
end 

c  Part  13:  top  water. 


start 
1  24; 
1  22: 

1   24  54  77: 
0.0  46.0 
4.0  46.0 
-106.0  -60.0 
nrb  2  0  0   2 
nrb  0  2  0 
nrb  0  0  1 
nrb  0  0  4 
mate  2 

end 


0.0  46.0 

0  0 

0  2  0 

0  0  1 

0  0  4 


c  Part  14:  side  water 


start 

1  22: 

1   3: 

1   24  54  77: 

4.0  46.0 

0.0  4.0 

-106.0  -60.0 

nrb  2  0  0 

nrb  0  0  1 

nrb  0  0  4 

mate  2 
end 


0.0  46.0 
2  0  0 
0  0  1 
0  0  4 


end 
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APPENDIX  B:  GRUNEISEN  GAMMA  APPROXIMATION 

The  seventh  order  approximation  for  the  Gruneisen  gamma  as  a 
function  of  the  specific  volume  v  developed  by  Gurtman.  Kirsch,  and 
Hastings  (1971)  is: 

y(v)  =  ao  +  ajv +  a2V2  +  •  •    -+a:v' 

where  a0  =  2.366.6324 

ai  =-22,669.420 
a2  =  91,259.368 
a3= -200, 175.85 
a4  =258,585.11 
a5  = -196, 872.84 
afi  =  81,850.023 
a7  =-14,342.530 
y  =  Gruneisen  gamma,  dimensionless 

v  =  specific  volume  in  i^f  • 

By  definition. 

P      ,      vo      , 
^pT"1^"1 

V0  1         r  ,        g  1 

giving  v  = r 7  tor  p  -  1  — -  =  77-. 

5        &  H  +  l      u+1        p        cm3      vo 

Substituting  l/(n+l)  for  v  into  the  Gurtman  equation  for  u  =  0  to  0.8  in 
.001  increments,  then  using  a  least  squares  linear  fit  with  the  same  intercept 
resulted  in  the  following  linear  equation  for  y(f.i): 

y(n)  =  Yo  +  ay  =  0.4934  +  1.3937u 
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Figure  41  is  a  plot  of  y(u)  from  the  Gurtman  equation  and  the  linear 
approximation. 
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Figure  41.    Comparison  of  y{\x)  from  Gurtman,  Kirsch,  and  Hastings 
equation  with  the  linearized  equation  for  y(n). 


62 


APPENDIX  C:  FORTRAN  PROGRAM 

This  appendix  contains  a  complete  listing  of  the  FORTRAN  main 
program  DTAPOPT  and  its  supporting  subroutines  DCOMPAR.  DPTGEN, 
and  DTAPWT,  written  to  optimize  tapered  charge  design.  A  separate 
subroutine  to  generate  pressure-time  histories  given  tapered  charge 
geometry  and  standoff  distance  is  required,  as  well  as  the  optimization 
package  ADS. 

PROGRAMMING  NOTE:  The  program  DTAPOPT  and  associated 
subroutines  were  written  in  double  precision  FORTRAN.  The  public  domain 
version  of  ADS  is  written  in  single  precision.  ADS  was  converted  to  double 
precision  for  use  with  DTAPOPT.  DTAPOPT  and  its  three  subroutines  can 
be  simply  converted  to  single  precision  by  removing  the  IMPLICIT  NONE 
and  DOUBLE  PRECISION  statements  from  the  codes. 

I.   PROGRAM  DTAPOPT 

c  DTAPOPT 

c         This  FORTRAN  program  combined  with  subroutines  DCOMPAR. 

c  DPTGEN,  and  DTAPWT.  is  designed  to  optimize  tapered  charge  design 

c  using  the  public  domain  optimization  package  ADS  (converted  to 

c  DOUBLE  PRECISION  by  the  author) 

c         --When  coupled  with  a  separate  subroutine  (not  included) 

c  to  calculate  the  pressure- time  history  of  a  three-segment  tapered 

c  charge . 

c         Up  to  ten  charge  segments  may  be  used  with  appropriate 

c  modifications  to  the  input  files. 

c_-_______________.  ----------------- 

c  PRECISION:  DOUBLE 

c  INPUT  FILES:  TAPIN.DAT   Initial  Design. 

c  PROFIN.DAT  Desired  P-T  History. 

c  INTERRACTIVE  INPUT:  Optimization  options. 

c  OUTPUT  FILE:  Named  in  TAPIN.DAT,  plot  data  for  computed 

c  time  vs  computed  pressure  and  interpolated  pressure. 

c  SCREEN  OUTPUT:  Starting  and  Optimized  Designs. 

c  REQUIRED  SUBROUTINES:  DTAPWT   Computes  charge  weight. 

c  DPTGEN   Linearly  interpolates  to  give  desired 
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AUTHOR : 

MOST  RECENT  UPDATE 

REFERENCES: 


pressure  at  each  time  output  by  P-T  history  generator. 

DCOMPAR  Computes  the  square  root  of  the 
average  of  the  square  of  the  difference  between  the  computed  and 
interpolated  desired  pressure  at  each  time. 

ADS     Package  of  subroutines  which  perform 
the  optimization. 

P-T  HISTORY  GENERATOR. 

William  Earl  Miller  II 

6/16/92 

(1)  Vanderplaats,  G.  N..  ADS  -  A  FORTRAN 
Program  for  Automated  Design  Synthesis,  Version  1.10,  program 
instructions,  Naval  Postgraduate  School.  Montery,  California, 
May.  1985. 

(2)  Vanderplaats,  G.  N.,  Numerical 
Optimization  Techniques  for  Engineering  Design:   With  Applica- 
tions, McGraw-Hill  Publishing  Company.  1984. 

TEST  P-T  GENERATOR:    The  subroutine  SDTAPER.  was  converted  to 
FORTRAN  by  the  author,  based  upon  the  PASCAL  program  TAPER  version 
7/27/89  by  F.  A.  Costanzo. 

DATA  INPUT  FILES- -Two  required: 

TAPIN.DAT  contains  12  lines,  one  value  per  line  which  are  read 
into  the  program.   The  integer  and  decimal  data  need  only  be  in 
a  format  suitable  for  list-directed  input  assignment  to  INTEGER 
and  DOUBLE  PRECISION  data  types  respectively.   The  character  data 
must  be  in  proper  form  for  a  DOS  file  name. 

READ  TO 
DESCRIPTION 
number  of  charge  segments 
length  of  charge  segment  (ft) 


LINE 
NUMBER 

1 

2 

3 

4 

5 

6 

7 

8 

9 
10 


diameter  of  charge  segment  (in) 


VARIABLE 
NTAP 
TAP(l) 
TAP(2) 
TAP(3) 
DIAM(l) 
DIAM(2) 
DIAM(3) 

'  DIAM(4) 
RANGE 


DATA 

TYPE 
integer 
decimal 
decimal 
decimal 
decimal 
decimal 
decimal 
decimal 
decimal 
decimal 


standoff  (ft) 

nominal  length  of  time  history  (ms)  TLEN 
used  by  P-T  generating  subroutine 

11  pressure  cutout  time  (ms)  TCO         decimal 

used  by  P-T  generating  subroutine 

12  name  of  output  data  file  to  be     NAMFIL      character 

created 
PROFIN.DAT  contains  NDES+2  lines,  one  integer  value  (NDES)  on  the 
first  line,  two  decimal  values  on  each  of  the  remaining  lines. 
Times  must  be  in  ascending  order  starting  with  0.0  and  extending 
to  a  time  greater  than  TLEN  entered  in  TAPIN.DAT  above 


LINE 
NUMBER 
1 


READ  TO 
VARIABLE 
NDES 


DESCRIPTION 
total  number  of  desired  pressure- 
time  history  data  pairs  minus  one 

2  first  time   'first  pressure      TDES(O)  PDES(O) 
(must  be  zero) 

3  second  time   second  pressure     TDES(l)  PDES(l) 


DATA 

TYPE 

integer 

decimal 

decimal 


NDES+2   NDES  time   NDES  pressure    TDES(NDES)  PDES(NDES)  decimal 
INTERRACTIVE  KEYBOARD  INPUT- -three  integers  required 
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These  integers  control  the  optimization  method  used  by  ADS. 

See  Ref  1  for  more  complete  instructions,  Ref  2  for' theory. 
Combinations  0  5  7.  8  5  7.  and  9  5  7  are  recommended.   Others  may 
work  well  in  a  given  design  space.   Using  the  initial  design,  three 
or  four  runs  with  different  combinations  should  result  in  a  most 
improved  design  which  can  then  be  used  as  the  initial  design  for 
further  optimization.   The  three  integers  are: 


ENTRY 
ORDER 
first 


READ  TO 
VARIABLE 
ISTRAT 


DATA 
TYPE 
integer 


second 


third 


DESCRIPTION 
optimization  strategy  used  by  ADS 
0  Go  directly  to  the  optimizer 

6  Sequential  Linear  Programming 

7  Method  of  Centers  (Design  must  be  feasible) 

8  Sequential  Quadratic  Programming 

9  Sequential  Convex  Programming 

optimizer  to  be  used  by  ADS  IOPT      integer 

4  Method  of  Feasible  Directions 

5  Modified  Method  of  Feasible  Directions 

one  dimensional  search  options  IONED     integer 

5  Golden  Section  Method 

6  Golden  Section  Method  Plus  Polynomial  Interpolation 

7  Bounded  Polynomial  Interpolation 

8  Unbounded  Polynomial  Interpolation 


OUTPUT  FILE:  Named  in  TAPIN.DAT.  consists  of 
three  columns.   Format  is  list-directed  from 
COLUMN 

DESCRIPTION 
times  output  from  P-T  GENERATOR 
P's  from  P-T  GEN'R  at  each  time 


NUMBER 
one 
two 
three 


int.  P's  for  each  time  from  DPTGEN 


NTIME+1  rows  of 
DP  variables. 

WRITTEN  FROM 
VARIABLE 
TIME(I).  I=0.NTIME 
PRESS(I).  I=0.NTIME 
PCOMPAR(I).  I=0.NTIME 


c  SCREEN  OUTPUT:   Outputs  initial  and  final  design  values  plut  the 

c  number  of  calls  to  the  P-T  Historv  Generator. 

c  INITIAL  AND  FINAL  VALUES 

c  DESCRIPTION 

c  lengths  of  tapered  charge  segments  (ft) 

c  diameters  of  tapered  charge        (in) 

c  charge  weight                     (lb) 

c  standoff                          (ft) 

c  sqrt  of  average  sq  diff  of  Pdesign-Pdesired 

c  OPTIMIZATION  EFFICIENCY 

c  number  of  calls  to  P-T  generator 


VARIABLE 
TAP(I) ,I=1.NTAP 
DIAM(I) , I=0.NTAP 
WEIGHT 
RANGE 
OBJ 

NCALLS 


VARIABLES 

NAME 


DATA    COMMON      ASSIGNED 
TYPE  BY 


A(I.J), 

1=1, NRA 
DF(21) 
DIAM(I), 

I=1,NTAP+1 
G(I), 

I=l,NCON 
I 
IC(I), 

I=1,NC0N 
IDG(I), 

I=1.NC0N 


DP 
,J=1,NC0LA 
DP 


ADS 


ADS 
DP     OPTDPA    TAPIN.DAT, MAIN 

DP       -  MAIN 

I       -  MAIN 

I       -  ADS 

I       -  MAIN 


USED     DESCRIPTION  ADS 
BY  arg 

? 

ADS  only    Constr  grads   Y 


ADS  only 
MAIN, 

PTGEN'R 
ADS 

MAIN 
ADS  only 

ADS 


Obj  grads  Y 

Charge  diams  N 

Constraints  Y 

Local  indexing  N 

Gradient  ID  Y 

Constraint  ID  Y 
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IGRAD 

INFO 

IONED 

IOPT 

IPRINT 

ISTPxAT 

IWK(I). 

I=1,NRIWK 
NAMFIL 
NCALLS 
NCOLA 
NCON 
NDES 
NDV 
NGT 
NRA 
NRIWK 
NRWK 
NTAP 
NT  I  ME 
OBJ 
PCOMPAR(I). 

1=0. NDES 
PDES(I) 

1=0, NDES 
PRESS(I) 

I=0.NTIME 
RANGE 
SUMSQ 
TAP(I) 

1=1, NTAP 
TCO 
TDES(I) . 

1=0. NDES 
TIME(I) 

I=0.NTIME 
TLEN 
VLB(I) , 

1=1, NDV 
VUB(I) . 

1=1. NDV 
WEIGHT 
WK(I) , 

1=1. NRWK 
X(I), 

1=1, NDV 


I 
I 
I 
I 
I 
I 
I 

CH 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
I 
DP 
DP 


DP 

DP 
DP 
DP 

DP 
DP 

DP 

DP 
DP 

DP 

DP 
DP 

DP 


FORDPTGEN 


OPTI 
OPTI 

FORDPTGEN 


DP   FORDPTGEN 


OPTDPA 

OPTDP 
FORDCOMPAR 
OPTDPA 

OPTDP 
FORDPTGEN 

OPTDPA 

OPTDP 


PASS 


MAIN 
MAIN, ADS 

user 

MAIN 

MAIN 

user 
ADS, MAIN 

TAPIN.DAT 
MAIN 
MAIN 
MAIN 
PROFIN.DAT 
MAIN 
ADS 
MAIN 
MAIN 
MAIN 
TAPIN.DAT 
PTGEN'R 

MAIN 
DPTGEN 

PROFIN.DAT 

PTGEN'R 


ADS 
MAIN. ADS 
ADS 
ADS 
ADS 
ADS 
ADS 

MAIN 

MAIN 

ADS 

ADS 
DPTGEN 

ADS 

ADS 

ADS 

ADS 

ADS 
PTGEN'R, MAIN 
DPTGEN . DCOMPAR 
MAIN  ADS 


Grad  Calc  Ctrl 
Prog  Flow  Ctrl 
One-D  Search 

Optimizer 
ADS  Print  Opts 
Strategy 
Work  Array 


DCOMPAR 


DPTGEN 


Output  File 
PTGEN'R  calls 
No.  A  columns 
No.  constraints 
Des'd  plot  sts-1 
No.  Design  vars 
Gradient  Ctrl 
A  rows 
IWK  Dimension 
WK  Dimension 
No.  Chg  Segs 
Pit  sts-1 
Obj  Func  Val 
Interp'd  Press 


PCOMPAR.MAIN 


Desired  Press  N 
Calc  Press   N 


TAPIN.DAT. MAIN   PTGEN'R       Standoff 

DCOMPAR        MAIN       Press  var'n 
TAPIN.DAT. MAIN  PTGEN'R, MAIN   Chg  Seg  Lth 


TAPIN.DAT 
PROFIN.DAT 


PTGEN'R 
DPTGEN 


Surf  CO  Time 
Des'd  PT  Time 


N 
N 
N 

N 

N 


PTGEN'R  DPTGEN. DCOMPAR. MAIN  Calc  time  N 


TAPIN.DAT 
MAIN 

MAIN 

DTAPWT 
ADS,  MAIN 

MAIN,  ADS 


PTGEN'R 
ADS 

ADS 

MAIN 
ADS 

ADS 


Lngth  PT  hist  N 
L  lim  on  DV  Y 

U  lira  on  DV  Y 


Charge  Weight 
Work  Array 

Des  Vas 


SUBROUTINES 
ADS ( INFO , I STRAT , IOPT , IONED , IPRINT , IGRAD . NDV , NCON . X . VLB . VUB , OBJ , 
G , I DG , NGT , I C , DF , A , NRA , NCOLA , WK , NRWK , IWK , NRIWK ) 
See  Ref  1  for  more  complete  instructions.  Ref  2  for  theory. 
Simply,  ADS  inputs  design  variables,  constraints,  and  the  objective 
function  for  an  initial  design,  then  modifies  that  design, 
requesting  the  corresponding  objective  and  constraint  values  from 
the  calling  program. 
ARGUMENT  VARIABLES 
A(NRA, NCOLA)  DP  Array  of  constraint  grads .   ADS  use  only  here. 
DF(NDV+1)     DP  Array  of  objective  gradients.  ADS  use  only  here. 
G(NCON)      DP  Array  of  constraints  for  current  design  in  X 
G(l&2):  25#<WEIGHT<125#;  G(3-5):  DIAM( 1)<DIAM(2 )<DIAM( 3)<DIAM(4) 
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IC(NGT)  I 
IDG(NCON)  I 
IGRAD  I 
INFO  I 

IONED  I 
IOPT  I 

IPRINT  I 
ISTRAT  I 
IWK(NRIWK)  I 
NCOLA  I 
NCON  I 
NDV  I 

NGT  I 

NRA  I 

NRIWK  I 
NRWK  I 
OBJ  DP 

VLB(NDV+1) 
VUB(NDV+1) 
WK 

X(NDV+1) 

X(I)=TAP(I) 

VLB(I)=1" 

VUB(I)=10' 

DCOMPAR 

Communicates  via 

Input  variables: 

Output  variable: 

DPTGEN 

Communicates  via 

Input  variables: 


DP 
DP 
DP 
DP 


Array  of  constraint  ID's.   NA  this  program 

Array  ID'ing  type  of  constraints:  0  for  nonlin  ineq 

=0  for  ADS  calculate  gradients  using  FD . 

ADS  flow  control  parameter. 

=5,6,7.8:   See  input  section. 

=4.5:  See  input  section. 

=0000  FOR  NO  ADS  PRINTOUT 

=0,6,7,8,9:   See  input  section 

Stores  ADS  I  vars .   Some  modifiable. 

Dimensioned  columns  of  A.  min  NDV+1 . 

Number  of  constraints  in  G. 

Number  of  design  variables  in  X. 

Returned  to  by  ADS  for  gradients.   0  this  program. 

Dimensioned  A  rows:   at  least  NDV+1. 

Est:  200+NDV+NCON+N+MAX(N,2*NDV) ,N=MAX( NDV. NCOLA) 

Est:  500+10*(NDV+NCON)+NCOLA*(NCOLA+3)+N*(N/2)+l 

Objective  function  value,  SUMSQ  provided  by  DCOMPAR 


Array  of  design  variable  lower  bounds,  indices  as 
Array  of  design  variable  upper  bounds,  indices  as 
Array  for  ADS  double  precision  variables. 
Arrav  of  design  vars.  Assigned  bv  input,  then  ADS, 
1=1,3/      X(4+I)=DIAM(I) .  1=0,3"     ;  X(8)=RANGE 
VLB(4+I)=0.75"  VLB(8)=5' 

VUB(4+I)=10"  VUB(8)=30' 


X 
X 


COMMONS 
NTIME: 
SUMSQ 

COMMONS 
NTIME: 
NDES; 


OPTI.  OPTDPA.  FORDPTGEN.  FORDCOMPAR 
PCOMPAR(I).  PRESS(I),  1=0. NTIME 


OPTI.  OPTDPA.  FORDPTGEN 
TIME(I),  1=0. NTIME 
TDES(I).  PDES(I).  1=0. NDES 


OPTDPA,  PASS 
I=0,NTAP.  TAP (I) 


I=1.NTAP 


Output  variables:  PCOMPAR(I).  1=1. NTIME 
DTAPWT 

Communicates  via  COMMONS:   OPTI. 

Input  variables:    NTAP :  DIAM(I) 

Output  variables:   WIEGHT 
PRESSURE -TIME  HISTORY  GENERATOR 

Communicates  via  COMMONS:   OPTI,  OPTDP .  OPTDPA 

Input  variables:   NTAP:  RANGE:  TLEN :  TCO ;  TAP(I).  1=1. NTAP:  DIAM(I) 
1=0, NTAP. 

Output  variables:   NTIME;  TIME(I).  PRESS(I).  1=0, NTIME 


PROGRAM  DTAPOPT 
IMPLICIT  NONE 

SPECIFICATIONS  FOR  SUB  VARIABLES 
COMMON /OPT I/NTAP , NTIME 
INTEGER     NTAP, NT I ME 
COMMON/OPTDP/    RANGE . TLEN , TCO 
DOUBLE  PRECISION  RANGE . TLEN . TCO 

COMMON/OPTDPA/    TAP.     DIAM,       TIME,         PRESS 
DOUBLE  PRECISION  TAP(10) , DIAM(0 : 10) , TIME(0 : 1000) . PRESS(0 : 1000) 
COMMON/PASS/     WEIGHT 
DOUBLE  PRECISION  WEIGHT 

COMMON/FORDPTGEN/NDES . TDES ,        PDES .        PCOMPAR 
INTEGER  NDES 

DOUBLE  PRECISION      TDES(0 : 1000) , PDES(0 : 1000) , PCOMPAR(0 : 1000) 
COMMON/FORDCOMPAR/SUMSQ 
DOUBLE  PRECISION   SUMSQ 


67 


SPECIFICATIONS  FOR  ADS  VARIABLES 
up  to  20  X's,  100  G's.  30  const  grads 
INTEGER  IWK(2000) . IDG(IOO) . IC(30) 

DOUBLE  PRECISION  X(21) , VLB(21) . VUB(21) . G(100) , DF(21) ,  A(21.30), 
c  WK(IOOOO) 

INTEGER  NRA . NCOLA , NRWK . NRIWK . IGRAD . NDV , NCON , ISTRAT , IPRINT . IOPT , 
c  IONED.  INFO,  NGT 
DOUBLE  PRECISION  OBJ 


SPECIFICATIONS  FOR  SUBROUTINES 


EXTERNAL  SDTAPER 
EXTERNAL  DTAPWT 
EXTERNAL  DPTGEN 
EXTERNAL  DCOMPAR 
EXTERNAL  ADS 

INTEGER  I.NCALLS 
CHARACTER*12  NAMFIL 


SPECIFICATIONS  FOR  LOCAL  VARIABLES 

BEGIN  EXECUTION 

ADS  ARRAY  DIMENSIONS 


NCALLS=0 

NRA=21 

NCOLA=30 

NRWK=10000 

NRIWK=2000 

ADS  PARAMETERS 
(no  user  provided  gradients,  5  constraints) 
(#  design  variables  determined  in  initial  design  section) 

IGRAD=0 

NCON=5 

INPUT  INITIAL  DESIGN.  OUTPUT  FILE, 
BOUND  ON  DESIGN  VARIABLES 

OPEN ( 88, FILE=' TAPIN.DAT' , STATUS=' OLD' ) 

READ(88.*)NTAP 


NDV=2*(NTAP+1) 

DO  99  I=1.NTAP 
READ(88.*)TAP(I) 
X(I)=TAP(I) 
VLB(I)=1.0D0/12.0D0 
VUB(I)=20.0D0 
99  CONTINUE 

DO  100  I=1.NTAP+1 
READ(88,*)DIAM(I-1) 
X(NTAP+I)=DIAM(I-1) 
VLB(NTAP+I)=0.75D0 
VUB(NTAP+I)=10.0D0 
100  CONTINUE 

READ(88.*)RANGE 
X(2*NTAP+2)=RANGE 
VLB(2*NTAP+2)=15.0D0 
VUB( 2*NTAP+2 )=30 . 0D0 

READ(88,*)TLEN 
READ(88,*)TCO 

READ(88,1000)NAMFIL 

CLOSE(88) 

OPEN(89,FILE=NAMFIL) 


1"  <  length  of  segment  <  20' 


0.75"  <  joint  diameter  <  10" 


5'  <  standoff  <  30' 

INPUT  LENGTH  OF  TIME  RECORD , TCO 
INPUT  DESIRED  PRESSURE  PROFILE 
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1000  FORMATCA12) 

OPEN C 88. FILE=' PROFIN.DAT' . STATUS=' OLD ' ) 

READ(88.*)NDES 

DO  101  I=0.NDES 

READ(88.*)TDES(I) .PDES(I) 
101  CONTINUE 
CLOSE(88) 

ID  FIVE  NONLINEAR  CONSTRAINTS 
IDG(1)=0 
IDG(2)=0 
IDG(3)=0 
IDG(4)=0 
IDG(5)=0 

INPUT 
(no  ADS  print,  interractive , optimizer  and  search  options) 
IPRINT=0000 

PRINT*.'  ENTER  ISTRAT . IOPT , IONED : 
READ* . I STRAT . IOPT , IONED 

OPTIMIZE 
(0  no  override.  -2  default  override) 
INFO— 2 
10  CALL  ADS(INFO. ISTRAT. IOPT. IONED. IPRINT . IGRAD ,NDV,NCON, X. VLB. 
C  VUB . OBJ . G . IDG . NGT . IC . DF . A . NRA . NCOLA . WK , NRWK . IWK , NRIWK) 

IFS  TO  CONTROL  FLOW  USING  ADS  INFO 


IF(INFO.EQ.l)  THEN 


DO  102  I=1.NTAP 
TAP(I)=X(I) 

102  CONTINUE 

DO  103  I=1.NTAP+1 
DIAM(I-1)=X(NTAP+I) 

103  CONTINUE 
RANGE=X(2*NTAP+2) 

CALL  SDTAPER 
NCALLS=NCALLS+1 

CALL  DTAPWT 

CALL  DPTGEN 

CALL  DCOMPAR 
OBJ=SUMSO 


EVALUATE  OBJECTIVE  AND  CONSTRAINTS 
ADS  VARIABLES  TO  SDTAPER  INPUT 


CALCULATE  PRESSURE  PROFILE 

CALCULATE  WEIGHT 

CORRESPONDING  DESIRED  PRESSURES 

EVALUATE  OBJECTIVE  FUNCTION 

PRINT  INITIAL  DESIGN  DATA  TO  SCREEN 


IF(NCALLS.EO.l)  THEN 

PRINT*.'       INITIAL  DESIGN' 
PRINT1001, (TAP(I) ,I=1.NTAP) 
PRINT1002 , ( DIAM( I ) . 1=0 , NTAP ) 
PRINT1003, WEIGHT 
PRINT1005, RANGE 
PRINT1004.0BJ 
ENDIF 


G(1)=(25.0D0-WEIGHT) 
G(2)=(WEIGHT-75.0D0) 

G(3)=X(4)-X(5) 
G(4)=X(5)-X(6) 
G(5)=X(6)-X(7) 


EVALUATE  CONSTRAINTS  (25//  <  W  <  125//) 
EVALUATE  CONSTRAINTS  (DKD2<D3<D4) 
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GO  TO  10 
ELSE 
IF(INF0.EQ.2)  THEN 
C  RESERVED  FOR  ADS  GRADIENTS 

GO  TO  10 
ELSE 

IF(INFO.EQ. -1)  THEN 
C  ADS  OVERRIDE  VALUES 

c  (for  no  scaling,  rel  FD  step,  min  | FD  step |, zero) 

IWK(2)=0 
WK(21)=0.001D0 
WK(22)=0.0001DO 
WK(26)=0.1D0 
WK(37)=1.0D-10 
GO  TO  10 
ELSE 
C  OPTIMIZATION  COMPLETE,  INFO=0 

IF(INFO.EQ.O)  GO  TO  20 
END  IF 
ENDIF 
ENDIF 
20  CONTINUE 
C  RECORD  RESULTANT  PRESSURE  PROFILES 

DO  104  I=0.NTIME 

WRITE(89.*)  TIME(I) .PRESS(I) .PCOMPAR(I) 
104  CONTINUE 
CLOSE(89) 
C  PRINT  DESIGN  DATA  TO  SCREEN 

PRINT* 

PRINT*. '       FINAL  DESIGN' 
PRINT1001. (TAP(I),I=1,NTAP) 

1001  FORMAT  ('         LENGTHS  =     ',3F10.3,'  FT') 
PRINT1002, (DIAM(I),I=0,NTAP) 

1002  FORMAT ('         DIAMETERS  =  ',4F10.2,'  IN') 
PRINT1003, WEIGHT 

1003  FORMAT ('         CHARGE  WEIGHT  =',F6.1,'  LB') 
PRINT1005, RANGE 

1005  FORMAT ('  RANGE  =     '.F10.1,'  FT') 
PRINT1004,OBJ 

1004  FORMAT ('         SQRT  OF  AVG  OF  (Pdesign-Pdesired) "2  =' .  F6.1, 
c'  PSI') 

PRINT* 
PRINT1006.NCALLS 

1006  FORMATC        CALLS  TO  P-T  HISTORY  GENERATOR  =',I4) 
END 


II.     SUBROUTINE  DCOMPAR 

c  SUBROUTINE  DCOMPAR 

c    This  subroutine,  for  use  with  DTAPOPT ,  computes  the  square  root  of 

c   of  the  average  sum  of  the  squares  of  pressure  differences. 

c----------------------------------- 

c   PRECISION:     DOUBLE 

c  AUTHOR:       William  Earl  Miller  II 

c   LAST  UPDATE:   6/17/92 

c   INTERFACE:     4  BUSES:   OPTI ,  OPTDPA,  FORDPTGEN .  FORDCOMPAR . 
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VARIABLES 
NAME 
I 

NT  I  ME 

PCOMPAR( I , 1=0 . NTIME) 
PRESSCI, 1=0, NTIME) 
SUMSO 


DESCRIPTION 

local  index 

calculated  t 

intp'd  desired  P 

calculated  P 

value  to  be  opt'd 


TYPE  COMMON  OR  LOCAL  I/O? 

I         local         na 

I  OPTI         I 

DP         FORDPTGEN       I 

DP         OPTDPA        I 

DP       FORDCOMPAR      0 

******************************************************************* 

SUBROUTINE  DCOMPAR 
IMPLICIT  NONE 

SPECIFICATIONS  FOR  GLOBAL  VARIABLES 
COMMON/OPT I /NTAP , NTIME 
INTEGER     NTAP, NTIME 

COMMON/OPTDPA/    TAP,     DIAM,       TIME.         PRESS 
DOUBLE  PRECISION  TAP( 10) . DIAM(0 : 10) , TIMEfO : 1000) . PRESS(0 : 1000) 
COMMON/FORDPTGEN/NDES , TDES ,         PDES ,         PCOMPAR 
INTEGER  NDES 

DOUBLE  PRECISION      TDES(0 : 1000) , PDES(0 : 1000) . PCOMPAR(0 : 1000) 
COMMON/FORDCOMPAR/SUMSO 
DOUBLE  PRECISION   SUMSQ 

SPECIFICATIONS  FOR  LOCAL  VARIABLE 
INTEGER  I 

BEGIN  EXECUTION 
SUMSQ=0 . 0D0 
DO  99  1=0, NTIME 

SUMSQ=SUMSO+( PCOMPAR ( I ) -PRESS ( I ) )*( PCOMPAR( I ) -PRESS ( I ) ) 
99  CONTINUE 

SUMSO=SQRT( SUMSO/ (NTIME+1 ) ) 

RETURN 

END 


III.  SUBROUTINE  DPTGEN 


c  SUBROUTINE  DPTGEN 

c       This  subroutine,  for  use  with  DTAPOPT.  interpolates  the  input 

c  desired  pressure- time  history  to  find  the  desired  pressure  at  each 

c  time  computed  by  the  P-T  generator. 


c 

PRECISION: 

DOUBLE 

c 

AUTHOR : 

William 

Earl  Mi 

Lller  II 

c 

LAST  UPDATE 

6/17/92 

c 

INTERFACE: 

3  BUSES: 

OPTI 

OPTDPA.  FORDPTGEN 

c 
c 

VARIABLES 

c 

NAME 

TYPE 

COMMON  OR  LOCAL 

I/O? 

DESCRIPTION 

c 

I 

I 

local 

na 

local  index 

c 

J 

I 

local 

na 

local  index 

c 

NDES 

I 

FORDPTGEN 

I 

No.  input  pts  -  1 

c 

NTIME 

I 

OPTI 

I 

No.  calc  pts  -  1 

c 

PCOMPAR(I) 

1=0, NTIME 

DP 

FORDPTGEN 

0 

Int  Desired  Press 

c 

PDES(I),  1= 

=0 

.NDES 

DP 

FORDPTGEN 

I 

Input  pressure 

c 

TDES(I),  1= 

-0 

,NDES 

DP 

FORDPTGEN 

I 

Input  time 

c 

TIME(I),  1= 

-0 

, NTIME 

DP 

OPTDPA 

I 

Computed  time 

c  ******************************************************************** 


SUBROUTINE  DPTGEN 
IMPLICIT  NONE 

COMMON/OPT I /NTAP , NTIME 


SPECIFICATIONS  FOR  GLOBAL  VARIABLES 
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13). AND. (ABSfTDES(J)) . GE . 1 . OD-13) )  THEN 
PROFILE  START  TIME  NOT  ZERO' 


100 


INTEGER     NTAP.NTIME 

COMMON/OPTDPA/    TAP,     DIAM.       TIME.         PRESS 

DOUBLE  PRECISION  TAP(IO) .DIAMCO : 10) . TIME(0 : 1000) , PRESS(0 : 1000) 

COMMON/FORDPTGEN/NDES . TDES .         PDES .         PCOMPAR 

INTEGER  NDES 

DOUBLE  PRECISION      TDES(0 : 1000) , PDES(0 : 1000) . PCOMPAR(0 : 1000) 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
INTEGER  I.J 

BEGIN  EXECUTION 
J=0 

IF((ABS(TIME(J)) .GE.1.0D 
PRINT* , ' INPUT  PRESSURE 
STOP 
ELSE 
TIME(J)=TDES(J) 
END  IF 
DO  100  I=0,NTIME 

IF(TIME(I).GE.TDES(J)) 

IF(TIME(I).LE.TDES(J+1))  THEN 

PCOMPAR(I)=PDES(J)+((TIME(I)-TDES(J))/(TDES(J+l)-TDES(J))) 
:  *(PDES(J+1)-PDES(J)) 

GO  TO  100 
ENDIF 
END  IF 
J=J+1 
IF(J.EQ.NDES+1)  THEN 

PRINT*. 'OUT  OF  INPUT  PLOT  POINTS' 
STOP 
ELSE 
GO  TO  1 
ENDIF 
CONTINUE 
RETURN 
END 


THEN 


IV.  SUBROUTINE  DTAPWT 


c 
c 

SUBROUTINE  DTAPWT 

This  subroutine, 

for  use  with  DTAPOPT ,  cal 

culates  the  weight 

c 
c 
c 

of  a  tapered  charge . 

PRECISION:     DOUBLE 

c 

AUTHOR:       William 

Earl  Miller  II 

c 

LAST  UPDATE:   6/17/92 

c 

INTERFACE:     3  BUSES 

:  OPTI,  OPTDPA,  PASS 

c 
c 

VARIABLES  FROM  COMMON 

STATEMENTS 

c 

NAME             TYPE 

COMMON  OR  LOCAL   I/O? 

DESCRIPTION 

c 

DIAM(I),I=0,NTAP  DP 

OPTDPA          I 

Chg  Diams  (in) 

c 

I                I 

local         na 

local  index 

c 

NTAP             I 

OPTI          I 

No.  Chg  Segs 

c 

PI               DP 

local         na 

Pi 

c 

RHO              DP 

local         na 

H20  dens  (lbm/ft"3) 

c 

SGRAV            DP 

local         na 

Sp  Grav  of  Chg 

c 

TAP(I) ,1=1, NTAP   DP 

OPTDPA         I 

Length  of  Seg  (ft) 

c 

VOL              DP 

local         na 

Charge  Vol 

c 

WEIGHT           DP 

PASS          na 

Charge  Wt  (lbm) 

c 

*************************** 

SUBROUTINE  DTAPWT 
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IMPLICIT  NONE 

SPECIFICATIONS  FOR  GLOBAL  VARIABLES 
COMMON/OPT I /NTAP , NT I ME 
INTEGER     NTAP.NTIME 

COMMON/OPTDPA/   TAP.     DIAM,       TIME,        PRESS 
DOUBLE  PRECISION  TAP(IO) ,DIAM(0 : 10) ,TIME(0 : 1000) , PRESS(0 : 1000) 
COMMON/PASS/     WEIGHT 
DOUBLE  PRECISION  WEIGHT 

SPECIFICATIONS  FOR  LOCAL  VARIABLES 
DOUBLE  PRECISION  PI ,RHO , SGRAV, VOL 
INTEGER  I 

BEGIN  EXECUTION 
PI=4 . 0D0*ATAN(1 . 0D0) 
RHO=  62.32D0 
SGRAV=  1.712D0 
VOL=0 . 0D0 
DO  101  1=1, NTAP 

VOL=VOL+PI*TAP( I )*(DIAM( I  - 1 )*DIAM( I -1)+DIAM( I )*DIAM( I ) 
c  +DIAM(I-1)*DIAM(I))/1728.0D0 
101  CONTINUE 

WE I GHT=VOL*SGRAV*RHO 

RETURN 

END 
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